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Abstract. We investigate spherically symmetric cosmological models in Einstein-aether the¬ 
ory with a tilted (non-comoving) perfect fluid source. We use a 1+3 frame formalism and 
adopt the comoving aether gauge to derive the evolution equations, which form a well-posed 
system of first order partial differential equations in two variables. We then introduce nor¬ 
malized variables. The formalism is particularly well-suited for numerical computations and 
the study of the qualitative properties of the models, which are also solutions of Horava 
gravity. We study the local stability of the equilibrium points of the resulting dynamical 
system corresponding to physically realistic inhomogeneous cosmological models and astro- 
physical objects with values for the parameters which are consistent with current constraints. 
In particular, we consider dust models in (/?—) normalized variables and derive a reduced 
(closed) evolution system and we obtain the general evolution equations for the spatially 
homogeneous Kantowski-Sachs models using appropriate bounded normalized variables. We 
then analyse these models, with special emphasis on the future asymptotic behaviour for 
different values of the parameters. Finally, we investigate static models for a mixture of a 
(necessarily non-tilted) perfect fluid with a barotropic equations of state and a scalar field. 
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1 Introduction 


Since the vacuum in quantum gravity may determine a preferred rest frame at the microscopic 
level, gravitational Lorentz violation has been studied within the framework of general rela¬ 
tivity (GR), where the background tensor field(s) breaking the symmetry must be dynamical 
[1]. Einstein-aether theory [2, 3] consists of GR coupled, at second derivative order, to a 
dynamical timelike unit vector field, the aether. In this effective field theory approach, the 
aether vector field and the metric tensor together determine the local spacetime structure. 

The aether spontaneously breaks Lorentz invariance by picking out a preferred frame 
at each point in spacetime while maintaining local rotational symmetry (breaking only the 
boost sector of the Lorentz symmetry). Since the aether is a unit vector, it is everywhere 
non-zero in any solution, including flat spacetime. A systematic construction of an Einstein- 
aether gravity theory with a Lorentz violating dynamical field that preserves locality and 
covariance in the presence of an additional ‘aether’ vector field has been presented [2-6]. 

In the infra-red limit of (extended) Horava gravity [7] [a candidate ultra-violet comple¬ 
tion in the consistent non-projectable extension of Horava-Lifschitz gravity], the aether vec¬ 
tor is assumed to be hypersurface-orthogonal; hence every hypersurface-orthogonal Einstein- 
aether solution is a Horava solution (most of the solutions studied). The relationship between 
Einstein-aether theory and Horava gravity is further clarified in [8], where it is shown how 
Horava gravity can formally be obtained from Einstein-aether theory in the limit that the 
twist coupling constant goes to infinity. 

Cosmological models in aether theories of gravity are currently of interest. The impact 
of Lorentz violation on the inflationary scenario has been explored [4-6] (also see the review 
[9]). 1 In particular, the primordial spectra of perturbations generated by inflation in the 
presence of a timelike Lorentz-violating vector field has been computed, and the amplitude 
of perturbation spectra were found to be modified which, in general, leads to a violation of 
the inflationary consistency relationship [10]. 

In particular, it is of importance to study inhomogeneous cosmologies, in both GR and 
alternative gravitational theories, partially motivated by current cosmological observations. 
Measurements of anisotropies of the cosmic microwave background (CMB) from experiments 
including the WMAP [11] and Planck [12] satellites, have provided strong support for the 
standard model of cosmology with dark energy (and specifically a cosmological constant, A). 
However the latest measurements are in tension with local measurements of the Hubble ex¬ 
pansion rate from supernovae la [13] and other cosmological observables which point towards 
a lower growth rate of large-scale structure (LSS) (which may be evidence for deviations from 
the standard ACDM cosmological model). The possible observation by the BICEP2 exper¬ 
iment [14] of B-rnode polarisation in the CMB in excess of the signal due to lensing would 
indicate the first detection of gravitational waves, perhaps generated during an inflationary 
era. In particular, there is a growing body of work on the imprints of gravitational waves on 
large-scale structure. 

In this paper we will study spherically symmetric Einstein-aether models. We shall 
study perfect fluid matter models in general, and various subcases in particular. In a com¬ 
panion paper [15] we study spherically symmetric Einstein-aether scalar field models with an 

1 We note that in scalar field models in which the dimensionless parameters of the models are not constant 
(e.g., depend on the scalar field), it was found that inflationary solutions are possible even in the absence of 
a scalar field potential [5]. 
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exponential self-interaction potential. Einstein-aether models with an exponential potential 
were recently studied in [16-18]. 

We shall use the 1+3 frame formalism [19-21] to write down the evolution equations for 
non-comoving perfect fluid spherically symmetric models and show they form a well-posed 
system of first order partial differential equations (PDEs) in two variables. We adopt the so- 
called comoving aether gauge (which implies a preferred foliation, the only remaining freedom 
is coordinate time and space reparameterization). We introduce normalized variables. The 
formalism is particularly well-suited for numerical and qualitative analysis [22], 

In particular, we derive the governing equations for an aether and a tilted perfect fluid 
assuming that the acceleration u is non-zero and introduce (so-called /3~) normalized variables 
(some of the technical details are relegated to the Appendix B). The evolution equations are 
presented in various different forms. We also rigorously derive the evolution equations when 
u = 0. We also consider the special subset U = v = 0 (where U is the normalized acceleration, 
v is the tilt, and we also assume the model parameters eg = 0 and c a / 0) and derive the final 
reduced phase space equations in normalized variables. We briefly review the Friedmann- 
Lemaitre-Robertson-Walker (FLRW) models in which the source must be of the form of a 
comoving perfect fluid (or vacuum) and the aether must be comoving. We study in detail a 
number of special cases of particular physical interest. 

We first consider dust models, which are of particular interest at late times. We in¬ 
vestigate a special dust model with U = 0 and v = 0 in normalized variables (assuming 
c a / 0) and derive a reduced (closed) evolution system. The FLRW models in this spe¬ 
cial dust model correspond to an equilibrium point. We are particularly interested in the 
future asymptotic behaviour of the models for different values of the parameters. We then 
consider the spatially homogeneous Kantowski-Sachs models using appropriate normalized 
variables (non- (5 —normalized variables which are bounded), and obtain the general evolution 
equations. A full global dynamical analysis of these models is possible. We then consider 
a special case and analyse the qualitative behaviour for physically reasonable values of the 
parameters at both early and late times. Finally, we consider static models for a mixture of 
a (necessarily non-tilted) perfect fluid with a barotropic equations of state and a scalar field, 
which are also of physical importance (although perhaps more from the astrophysical point 
of view than from the cosmological one). A brief discussion of the physical conclusions is 
presented at the end. 

1.1 The models 

The evolution equations follow from the field equations (FE) derived from the Einstein-aether 
action [2, 3]. In an Einstein-aether model there will be additional terms in the FE which 
include (see the technical details in the Section 2.2): 

• The effects on the geometry from the anisotropy and inhomogeneities (e.g., the curva¬ 
ture) of the spherically symmetric models under consideration. 

• The Einstein FE are generalised by the contribution of an additional stress tensor, T®?, 
for the aether field which depends on the dimensionless parameters of the aether model 
(e.g., “the Cj”). In GR, all of the Cj = 0. To study the effects of matter, we could 
perhaps assume the corresponding GR values (or close to them) in the first instance. 

• When the phenomenology of theories with a preferred frame is studied, it is generally 
assumed that this frame coincides, at least roughly, with the cosmological rest frame 


- 3 - 


defined by the Hubble expansion of the universe. In particular, in an isotropic and 
spatially homogeneous Friedmann universe the aether field will be aligned with the 
(natural preferred CMB rest frame) cosmic frame and is thus related to the expansion 
rate of the universe. In principle, the preferred frame determined by the aether can be 
different from (i.e., tilted with respect to) the CMB rest frame in spherically symmetric 
models. This adds additional terms to the aether stress tensor Tj|, which can be 
characterized by a hyperbolic tilt angle, v(t), measuring the boost of the aether relative 
to the (perfect fluid) CMB rest frame [4, 5]. The tilt is expected to decay to the future 
in anisotropic but spatially homogeneous models [23]. 

1.2 Spherical symmetry 

All spherically symmetric aether fields are hypersurface orthogonal and, hence, all spherically 
symmetric solutions of aether theory will also be solutions of the IR limit of Horava gravity. 
The converse is not true in general, but it does hold in spherical symmetry for solutions with 
a regular center [ 8 ] . 

The Ci are dimensionless constants in the model. When spherical symmetry is imposed 
the aether is hypersurface orthogonal, and so it has vanishing twist. Thus it is possible to 
set C 4 to zero without loss of generality [1]. After the parameter redefinition to eliminate 
C 4 , one is left with a 3- dimensional parameter space. The c* contribute to the effective 
Newtonian gravitational constant G; so a renormalization of the parameters in the model can 
be then used to set 87 tG = 1 (i.e., another condition on the q can effectively be specified). 
The remaining parameters in the model can be characterized by two non-trivial constant 
parameters. The other constraints imposed on the q have been summarized in [1] (e.g., see 
equations 43-46 in [24]; also see Appendix A). In GR q = 0. We shall study the qualitative 
properties of models with values for the non-GR parameters which are consistent with current 
constraints. 

Some of the models studied in this paper involve a static metric coupled to a stationary 
aether. This situation will be referred to here as “stationary spherical symmetry”. This 
case will be treated separately later. An important special case occurs when the aether is 
parallel to the Killing vector. We refer to this special case as a “static aether”. A spherically 
symmetric static vacuum solution is known explicitly [25]. 

1.3 Stars and black holes 

Spherically symmetric static and stationary solutions are physical important. Unlike GR, 
Einstein-aether theory has a spherically symmetric mode, corresponding to radial tilting 
of the aether. The time-independent spherically symmetric solutions and black holes were 
studied in [25] and [26], respectively, and surveyed in [1], and recently revisited for a more 
viable coupling parameter ci in [25] and [24]. In general, within this same parameter space, 
the dynamics of the cosmological scale factor and perturbations differ little from GR, and 
non-rotating neutron star and black hole solutions are quite close to those of GR. A thorough 
examination of the fully nonlinear solutions has not been carried out to date. A fully nonlinear 
energy positivity has, however, been established for spherically symmetric solutions at a 
moment of time symmetry [27]. 

Let us discuss this in more detail. There is a three-parameter family of spherically 
symmetric static vacuum solutions [25]. In the Einstein-aether theory the aether vector 
and its derivative provide two additional degrees of freedom at each point. If asymptotic 
flatness is imposed and the mass is fixed, there remains a one-parameter family (i.e., imposing 
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asymptotic flatness reduces this to a two parameter family [28] ), whereas GR has the unique 
Schwarzschild solution (Birkhoff’s theorem). In GR asymptotic flatness is a consequence of 
the vacuum field equations without any tuning of initial data, so the one-parameter family 
of local (Schwarzschild) solutions is automatically asymptotically flat. The radial tilt of the 
aether provides another local degree of freedom in aether theory, so spherical solutions need 
not be time-independent (even when restricting to stationary spherically symmetric aether 
theory). Not only are spherical solutions not generally static, but even if we restrict to 
static, spherical solutions, they are not necessarily asymptotically flat. It was shown in [3] 
that the Reissner-Nordstrom metric in a spherically symmetric static gauge with fixed norm 
is a solution, although this is not the only solution in that special case [25]. 

Requiring that the aether be aligned with the timelike Killing field restricts the static 
aether solution to one parameter (the single parameter C 14 , essentially the total mass [25]). 
Thus the solution outside a static star is the unique vacuum solution for a given mass in the 
static aether case [25], and is asymptotically flat. In [29] it was found that this static “worm- 
hole” aether solution is generally stable to linear perturbations under the same conditions as 
for flat spacetime. In the pure GR limit (ci = 0), we have just the Schwarzschild solution. For 
small values of r, the solutions can behave quite differently from the Schwarzschild solution. 
More recently, an analytic static spherically symmetric vacuum solution in the Einstein-aether 
theory was presented (demonstrated numerically) by use of the Euler-Lagrange equations [30] . 

Unlike the singular wormhole, the static solutions have a regular origin [25]. It is known 
that pure aether stars do not exist; i.e., there are no asymptotically flat self-gravitating aether 
solutions with a regular origin [25]. It has been shown that in the presence of a perfect fluid, 
regular asymptotically flat star solutions exist and are parameterized (for a given equation 
of state) by the central pressure (see also [31]). 

For black holes the aether cannot be aligned with the Killing vector, since the latter 
is not timelike on and inside the horizon. Instead, the aether is at rest at spatial infinity 
and flows inward at finite radii. The condition of regularity (at the spin-0 horizon) selects a 
unique solution from the one-parameter family of spherical stationary solutions for a given 
mass [25, 26]. Such black holes are rather close to Schwarzschild outside the horizon for a 
wide range of couplings. Inside the horizon the solutions differ more (but typically no more 
than a few percent), and like the Schwarzschild solution they contain a spacelike singularity. 

More recently, static spherically symmetric, asymptotically flat, regular (non-rotating) 
black-hole solutions in Einstein-aether theory have been studied (numerically) [24], gener¬ 
alizing previous results. It has been found that spherical black-hole solutions formed by 
gravitational collapse exist for all viable parameter values of the theory and a notion of 
black hole thus persists. Indeed, static spherically symmetric solutions in Lorentz-violating 
theories, in which the causal structure of gravity is greatly modified, still possess a special 
hypersurface, called a “Universal horizon”, that acts as a genuine absolute causal boundary 
because it traps all excitations, even those which could be traveling at arbitrarily high prop¬ 
agation speeds [24]. The Universal horizon satisfies a first law of black-hole mechanics [32], 
and evidence has been found that Hawking radiation is associated with the Universal horizon 
[33, 34], 

Finally, it would be of interest to determine the structure of rotating solutions; rapidly 
rotating black holes, unlike the non-rotating ones, might turn out to be very different from 
the Kerr metrics of GR. 
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2 Spherically symmetric Einstein-aether Models 

We shall use the 1+3 frame formalism [19, 20] to write down the evolution equations for 
spherically symmetric models as a well-posed system of first order PDEs in two variables. 
The formalism is particularly well-suited for studying perfect fluid spherically symmetric 
models [21], and especially for numerical and qualitative analysis [22]. We follow a similar 
approach to that in the resource paper [35] (wherein all relevant quantitites are explicitly 
defined). 

2.1 Restrictions on the kinematic and auxiliary variables: 

The metric is: 


ds 2 = —N 2 dt 2 + (ei l )~ 2 dx 2 + (e 2 2 )~ 2 (dd 2 + sin 2 ddp 2 ). (2.1) 

The Killing vector fields (KVF) are given by [36]: 

dp, cos ip d$ — sin ip cot id dp, sin p d$ + cos <p cot d dp. (2-2) 

The frame vectors in coordinate form are: 


e 0 = N l d t , ei = e^d x , e 2 = e 2 2 d$, e 3 = e 3 3 dp, 

where e 3 3 = e 2 2 /sinh. N, ei 1 and e 2 2 are functions of t and x. 

This leads to the following restrictions on the kinematic variables: 


(2.3) 


a a/ 3 = diag(—2<r + , <t+, o+), u) a/3 = 0, u a = (hi, 0,0), (2.4) 


where 

u Q = U^V/3U a -, 

u\ = ei In IV; 

on the spatial commutation functions: 


/ 0 0 n 13 \ 

a a = {a i,a 2 ,0), n a/3 = 0 0 0, 

\n 13 0 0 / 


where 


a\ = ei lne 2 2 , 

and on the matter components: 


= rt 13 = —-e 2 2 cot d] 


q a = (hi, 0,0), 7 x a p = diag(-27r + ,7r + ,7r + ). 


(2.5) 

( 2 . 6 ) 


(2.7) 


( 2 . 8 ) 


(2.9) 


The frame rotation Q a /3 is also zero. 

Furthermore, rii 3 only appears in the equations together with e 2 ni 3 in the form of the 
Gauss curvature of the spheres 


2 K \= 2(e 2 - 2n 13 )n 13 , 


( 2 . 10 ) 
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which simplifies to 


2 K = (e 2 2 ) 2 . (2.11) 

Thus the dependence on d is hidden in the equations. We will also use 2 K in place of e 2 2 . 
To simplify notation, we will write 


2 K, iii , «i 


as 


K, u , a. 


To summarize, the essential variables are 


IV, ei 1 , K, 8, cr+, a, ft, /j, gi, p, vr + , (2.12) 

where IV is the lapse function, ei 1 is the non null component of the frame vector ei, K is 
the Gauss curvature of the spheres, 8 is the (volume) rate of expansion scalar, is related 
to the magnitude of the rate of shear tensor (a measure of the anisotropies present in the 
model), a is the radial component of the object (spatial commutation function) a a , ft is 
the acceleration, p denotes the total energy density scalar, q\ is a component of the total 
energy current density vector, p is the total isotropic pressure scalar and tt + is related to the 
magnitude of the total anisotropic pressure tensor [19, 20]. 

In the case of spherical symmetry in Einstein-aether theory one must be careful in 
choosing the gauge (an additional gauge condition). Normally, in GR, spherically symmetric 
coordinates are chosen so that the metric is simplified (e.g., a choice for N) or so that the 
fluid is comoving. Here we chose the aether vector field to be aligned with the timelike frame 
vector eo (the comoving aether gauge, and hence in general N (t, x) cannot be simplified 
any further). This may make comparisons with GR difficult in some special cases. Our 
formulation is perhaps better suited for fluids/matter and cosmology, although the static 
case is not necessarily aligned (see later). 

We note that the tilt is defined relative to matter; one important question is to investi¬ 
gate whether this tilt decays to the future. 


2.2 Einstein-aether theory 

The action for Einstein-aether theory is the most general generally covariant functional of 
the spacetime metric g a b and aether field u a involving no more than two derivatives (not 
including total derivatives) [1, 27]. The action is [1, 37]: 


S 


1 


d^x^f—g 


R - K ab cd V a u c V b u d + A ( u c u c + 1) + C r 


(2.13) 


where 

K ab c d = Cig ab g c d + C 2^c^d + c 3 + C/pu a u b g cd . (2-14) 

The action (2.13) contains an Einstein-Hilbert term for the metric, a kinetic term for the 
aether with four dimensionless coefficients q, and A is a Lagrange multiplier enforcing the 
time-like constraint on the aether. 2 The convention used in this paper for the metric 

2 We set the vector norm to unity in order to obtain a unit time-like aether. Comparing with [27, 39], the 
tensor K!^ n was rescaled by a factor of 2, ca was taken with the opposite sign, and J a b was redefined taking 
the opposite sign (i.e., the constant c’s here and A have been rescaled by a factor of 2). 


- 7 - 





signature is (—b++) and the units are chosen so that the speed of light defined by the 
metric g ab is unity and k 2 = 8itG = 1. The field equations from varying (2.13) with respect 
to g ab , u a , and A are given, respectively, by [39]: 

Gab = T T ah OT (2.15) 

\u h = X7 a J a b + c 4 UaV 6 u a (2.16) 

u a u a = -1. (2.17) 

Here G a b is the Einstein tensor of the metric g a b ■ T^ b OT is the total energy momentum tensor, 
Tab° J = T^b + T^ at , where T^ at is the total contribution from all matter sources. We shall 
omit T^ at for the moment (and add in later for perfect fluid and scalar field sources), and 
so we begin with the vacuum case {C m = 0) first with a non-trivial aether stress-energy T^ b 
(which we will refer to as “pure” Einstein-aether theory which is a theory of the spacetime 
metric g a b and a vector field (the “aether”) u a ). 

The quantities J a b , u a and the aether stress-energy T^ b are given by 


r m = -K ab mn v b u n 

(2.18a) 

u a — U V b'Ua 

(2.18b) 

T^b = 2ci(V a u c V 6 u c - V c u a V c u b ) 

— 2[V c (tt(' a e/ 6)) + V c (u <7(ab)) — ^ c{ u (aJ b ) )] 2c4'U a 'Ufe+ 

T 2A U a U b T gab^u 

(2.18c) 

where 

C u = -K ab cd V a u c V b u d , 

(2.19) 

is the Einstein-aether Lagrangian [40]. 

Taking the contraction of (2.16) with u b and with the induced metric h bc 

:= g bc + u b u c 


we obtain the equations 


A = —u b X7 a J b — c\ii a u a , (2.20a) 

0 = h bc V a J b + c 4 h bc u a V b u a . (2.20b) 

We shall use the equation (2.20a) as a definition for the Lagrange multiplier, whereas the 
second equation (2.20b) leads to a set of restrictions that the aether vector must satisfy. 

The Einstein FE, Jacobi identities and contracted Bianchi identities gives a system of 
partial differential equations on the frame and commutator functions, while (2.18c) defines 
the components of the energy momentum tensor and (2.20b) gives one extra equation for 
the aether. We choose a gauge in which the aether is aligned with eg, the comoving aether 
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temporal gauge (all that then remains is the time and space reparameterization freedom): 


. 3 


eolei 1 ) = -§(0 - 6(j + )ei 1 , 

e o (K) = -l(0 + 3a + )K, 

e o (0) — ei(it) = -§0 2 - 6a+ + (u — 2 a)u — ^(p + 3 p), 

e o(<7+) - ^(-u ~ a ) = ~ 6a + ~^( a + u)u~ ^I< + vr+, 

e 0 (a) = -§(0 + 3 a + )(a + u)- \q\, 

e 0 (/i) + ei(gi) = -0(p + p) + 2(o - u)gi - 6 ct + 7t + , 

eofe) + ei(p) = -§(20 - 3o+)gi - 2(3a - u - ei)7r + - u(p + p), 

Constraints: 


(2.21a) 

(2.21b) 

(2.21c) 

(2.21d) 

(2.21e) 

(2.21f) 

(2.21g) 


ei (In IV) = u, (2.22a) 

ei(ln/C) = 2a, (2.22b) 

H = 3H 2 -3crl + K-3a 2 + 2e 1 (a), (2.22c) 

gi = —6a<r + + §ei(0 + 3(j_)_), (2.22d) 

where 1 (p, p, q \, 7r + ) = (p s , p^, q m , 7r^) can be computed from (2.18c): 

Pas = (ci - c 4 )(2ei - 4a + 'd)u - §(ci + 3c 2 + c 3 )0 2 - 6(ci + c 3 )cr 2 , (2.23a) 

Pas = §(ci + 3c 2 + c 3 )(2e 0 + 0)9 - 6(ci + c 3 )ct 2 + §(c 4 - c 4 )u 2 , (2.23b) 

fe = -|(ci - c 4 )(0 + 3fj+)u - 2(ci - c 4 )e 0 (tt), (2.23c) 

TTae = §(d - c 4 )u 2 + 2(c 4 + c 3 )(e 0 + 0)cr + . (2.23d) 


The aether equation (2.20b) becomes (and is true regardless of whether u is zero or not) 


(ci - c 4 )e 0 (w) = — §(ci - c 4 )(0 + 3a + )u + 6(c 4 + c 3 )a<r + 

+ §(ci + 3c 2 + c 3 )ei(0) - 2(ci + c 3 )ei(a + ). (2.24) 


To simplify these expressions it is convenient to make a reparameterization of the aether 
parameters, analogous to the one given in [8]: 


Q? c 2 T (ci T c 3 )/3, Cfj Cl T c 3 , c^ c\ c 3 , Cq, c 4 c 4 , 


where the new parameters correspond to terms in the Lagrangian relating to expansion, 
shear, acceleration and twist of the aether. Since the spherically symmetric models are 
hypersurface orthogonal the aether held has vanishing twist and is therefore independent of 
the twist parameter c u (the coupling c 4 — c 3 does not occur in the held equations (only c 4 + c 3 
does) [8]; this is equivalent to being able to set c 4 = 0 [1]). 

A second condition on the Cj can effectively be specihed by a renormalization of the 
Newtonian gravitational constant G. The remaining parameters in the model can therefore 

3 Note that some degenerate cases, including the static case below, may not be easily included in this 
approach. 

4 Here, for example, /r = fi tot is the total energy density. 
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be characterized by two non-trivial constant parameters. The other constraints imposed on 
the Ci have been summarized in [1], 

It may be useful later to define c 2 = 1 — 2c a <1. In particular, some special cases of 
interest are (see the Appendix A): case A: c a = \(1 — c 2 ) > 0, c a = — c a < 0 ,cg = 0: 
case B(ii): c a = ^(1 — c 2 ) > 0,c a = —^(1 — c 2 ),cg = 0: case C: c a = ^(1 — (?) > 0, eg = 
-g(l - (?) < 0, [c a = 0] 

The Lagrangian (2.19) becomes 


C u = ~ ( c a ii 2 + cg9 2 + 6c CT cr 2 ) (2.25) 

and the aether energy components (p, p, q\ , 7r_|_) = (p s , p s , q m , it m ) become 

Pie = — cq6 2 — 6 c a cr+ — c a (u + 2ei — 4 a)u, (2.26a) 

Pae = -\c a ii 2 - 6c ct ct 2 + c g (2e 0 + 0)9 , (2.26b) 

fe = | c a {? + 2>a + )u + 2 c a e 0 (u), (2.26c) 

^ =-§c a u 2 + 2c a (eo + 0)a + . (2.26d) 

and the aether equation (2.20b) reads 

c a e 0 (u) = -§ c a (0 + 3 a + )u - 6c a aa + - ei(c e 0 - 2c a a + ). (2.27) 


Combining all of the above equations, and assuming u ? 0 (the special case u = 0 will 
be dealt with later), we obtain 


eo(ei 1 ) = -\{6 - 6cr+)ei 1 , 

e o (K) = ~l(0 + 3a + )K, 

(3cg + 2c a ) e\(9) 2 ■ (a , Q \ 

eo( " ) ^ 3c a (2c. - 1) = + 3 *+>- 

(c a + l)ei(ft) 1 2 (c a + 1) v? 2 (c a + 1) ail 6 (2c ff - 1) a\ 
eo(v) - tt. -——— = + ~77 -—rr- 77 ;-——-r 


(3 Cg + 1) 


( 3 C 0 + 1) (3 Cg + 1) 


(3 eg + 1) ’ 


_ /_ n (Ca + l)ei(u) , 1 2 , ( 3 C 0 + 1 )0 2 (1 — 2c a ) au 

eo(u+)- TTTTT, -IT - = ~^ a + + 2 a + + ttttt;-tv + vttt;-77—I- 


3( 2 c ct - 1 ) 


18 (2c a - 1) 3 (2Co- - 1) 


+ 


(5c a + 2) u 2 


+ 


K 


6 ( 2 c ff - 1 ) 2 ( 2 c ct - 1 ) 2 ( 2 c ct - 1 ) 

/ x . (3c@ + 2 c ct ) ei(0) _ 1 .x //3 Q x 

e o(a) H-————-— — 3(0 + u)(9 + 3a + ). 


Constraints: 


(2.28a) 

(2.28b) 

(2.28c) 

(2.28d) 

(2.28e) 

(2.28f) 


ei(ln N) = it, 
ei(ln K) = 2 a, 

Q /> <}/2 

ei(o) + c a ei(ii) = (3 eg + 1) 9 2 - - (2c a - 1) a\ - — + — + 2 c a au - -y- 

t \ (3ce + 1) ei(0) 

ei(<T+)-----= 3a<7 + . 


6c^ — 3 


Commutator: 


[e 0 , ei] = iie 0 -\(9- 6cT + ) e i- 


(2.29a) 

(2.29b) 

(2.29c) 

(2.29d) 

(2.30) 
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Integrability conditions. In the Einstein-aether analysis, one of the “field equations” 
is the spatial projection (with the induced metric h bc ) of the equation obtained by the con¬ 
traction of the velocity variation of the action. In many cases this equation does not involve 
the appropriate time derivatives, and hence this equation is not an evolution equation, but 
rather it is a constraint [41]. In the spherically symmetric case here it can be shown that the 
constraint is conserved and is compatible with all of the other (evolution) equations. 


3 Aether and a tilted perfect fluid 


The energy momentum-tensor for the matter field is 


TUSi = -2 


SC 


HI J- r n , 


(3.1) 


with p to be specified. In general, the 4-velocity vector u of the perfect fluid is not aligned 
with the vector eo of a chosen temporal gauge. In spherically symmetric models, u is allowed 
to be of the form 

u = T(e 0 + uei), T = (1 - u 2 )-i, (3.2) 

where v is the tilt parameter. We choose a linear equation of state for the perfect fluid: 


P=( 7 - 1 )£» 


(3.3) 


where 7 is a constant satisfying 1 < 7 < 2. 


Then we obtain for the tilted fluid: 


P 

P 

<7i 


71+ 


G a 


1 — v 


:P 


(7 — 1)(1 — v 2 ) + | 7 u 2 


1 — V 7 


p 


IP 

1 - 2 V 

1 — v z 

1 7/1 


3 1 — u 2 


(3.4a) 

(3.4b) 

(3.4c) 

(3.4d) 


where G± = 1 ± (7 — l)v 2 . Thus (the total) p, p, q\ and ir + are given in terms of fi and v. 
These are then substituted into the evolution and constraint equations. 

Assuming ii / 0 (and under the general conditions that 2co- — 1/0, 3 cq — 1 / 0) we 
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obtain: 


eoCei 1 ) = -\{9 - 6 cr + )ei 1 , 
e 0 (K) = —1(6* + 3a + )K, 


e 0 (u) - 


- + 3(7+) + 


eo(0) - 
+ 


(3 c g + 2 c a ) e x (6) 2 . /n , 0 _ N , ^c a \xv 

3 c a ( 2 c a - 1 ) 

(c a + l)ei(tt) _ 1 2 (c a + l)u 2 2 (c a + 1 ) ait 6 (2c ct - 1 ) ct 2 


(3 c e + 1) 

A ((7 ~ 2)n 2 -37 + 2) 
2 (3c0 + 1) (1 - v 2 ) 


= 4r+ , 

(3c0 + 1) 


c a ( 2 c CT - 1 ) (1 - o 2 ) ’ 

+ 


(3c 0 + 1) 


(3ce + 1) 


\ (c a + l)ei(n) , 1_ 2 , (3c 0 + l)6» 2 | (l-2c a )an i 

e o(cr+;- 7777 - 77— — —V(T+ + -cr + + —— - — + —— - — + 

3 (2co- — 1) 2 T 18 (2 c ct — 1) 3(200--!) 


+ 


(5c a + 2) u 2 


+ 


K 


+ 


A ((7 +1> 2 - 1 ) 


6 (2co- — 1) 2(2 C(T -1) 2 (2 c ct — 1) 6 {2c a - 1) (1 - v 2 ) ’ 

7 jjLV 


{2>c g + 2c a )e l (9) _ l 

e o(a) H-3(2^ _ T)- _ _ 3( a + u )(® + 3(J +) _ 

ei (A) v (2 — j) 7Aei(o) 2~/fiv 2 cr + 7 (v 2 - 3 ) fi6 7 fiva 

e o{9) -77-77- =-77-1-—-2- 


2 (2Ccr — 1) (1 — v 2 ) 


G _ 


G_ 


G_ 


3G- 


G _ 


_ e i (A) (l ~ ^ 2 ) 2 (7-1) _ v (2 — 7 ) ei (u) = 2n (l - v 2 ) a + | 


7AG- G_ 

n (l — o 2 ) (37 — 4) 6 2v 2 (l — v 2 ) ( 7 — 1)0 


G- 


3 G_ 


+ 




+ (1 — n 2 )n. 


(3.5a) 

(3.5b) 

(3.5c) 


(3.5d) 


(3.5e) 

(3.5f) 

(3.5g) 


(3.5h) 


Constraints: 


ei(lnlV) = u, 
ei(ln K ) = 2a, 

ei(o) + c a ei(ti) = 


( 2 *- 1 ) 4 + 


A' 3a 2 c a u 2 

-- + — + 2c a au~ — 

ei (<r + ) - (3 ?,„ +1)e . 1 A ) = 3a<r+ + 


7/70 


(3.6a) 

(3.6b) 


(3.6c) 

(3.6d) 


3(2cq- — 1) ' 2(2c CT -l)(l-o 2 )' 

3.1 Well-posedness 

We now show that the system of evolution equations plus restrictions for the state vector 

, T 


[ei\ K,u,6,a + ,a, fx,v\ 
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is: 


is well-posed for 7 > 1 . The coefficient matrix for the spatial derivative terms (for [ft, 9 , a + ,a, fi, v] A 
5 

(3^+2^) q q o 0 

3cq(2e CT -l) 

0 0 0 0 0 


0 

Cg + 1 


\ 


V 


(3c fl +l) 

Cg-h 1 

3(2c ct -1) 

0 

0 

0 


0 

3c/)+2c a 

3(2c a -l) 

0 

0 


0 0 
0 0 

0 0 
0 0 


0 

0 

v (2 - 7) 


0 

0 

tA 


(7 — l)v 2 — 1 (7 — l)v 2 — 1 

(1 -n 2 ) 2 (7- 1) v(2- 7 ) 


Its eigenvalues are 

0 , 0, =t 1 


((7 — l)v 2 — 1) 7 ft, (7 — l)w 2 — 1 J 
(c a + 1) ( 3 CQ + 2 c a ) (2 - 7)1; ± V7 - 1(1 - v 2 ) 


3 c a (3c0 + 1 ) (2c ct — 1) ’ 
with corresponding eigenvectors (for example) 


G_ 


/ 


1 


\ 


(Cn + l)(6Co-—3) 
Ca(3cfl + l)(3ce+2c g -) 
(Cq-|-l)(3cg + l) 
Ca(3c9+2c 0 -)(6c <r —3) 

1 

0 


V 
/ 0 \ 
0 
0 
1 
0 

Vo 


0 

0 

0 

0 

0 

IE 

( 1—1 

1 


/ 


1 


\ 


(c n + l)(6c CT —3) 
Cq(3Cfl + l)(3Cfl+2c 0 -) 
(Cq-t-l)(3cg + l) 

Ca (3cg+2c a ) (6 Cq -3) 

1 

0 


V 


0 

0 

0 

0 

0 

IE 

(1 — 1 

1 



/°\ 


0 


1 

5 

0 


0 


w 


( 3 . 7 ) 


( 3 . 8 ) 


( 3 . 9 ) 


The matrix is diagonalizable (all its eigenvalues are real) for 7 > 1 and JV^e+iXel: 2 ^) — 
0 , with c s = — 1 being the speed of sound in the perfect fluid. The system ( 3 . 5 a)-( 3 . 5 h) 

is thus well-posed if the above conditions are fulfilled, otherwise the system is elliptic and 
not well-posed. 


3.2 Normalized variables 

We introduce the normalized variables (for u 7^ 0 using the / 3 -normalization): 

{A f- 1 ,E 1 1 ,Q,E,A,u] = {N-\ell,a + ,a,u}//3 

{12, n k ,JC,S + } = {Cl,-± 3 R, K, 3 S+} /( 3 / 3 2 ), 
where f 3 = |(0 + 3 cr + ) and 

3 R = 4eia — 6a 2 + 2 K, 3 S + = — |eia + ^K. 

5 Strictly speaking, we should also include the factor ei 1 in the matrix, but the result on well-posedness is 
the same. 
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By definition Q + £ = 1. In the above we assume that /3 A 0. In general the variables are 
unbounded. However, physically 17 > 0, and if the expansion and the shear are both positive, 
then 0 < <2 < 1. 

We also introduce the normalized differential operators, 


d ° := 7 


; where a = 0,1. 


Moreover, we define q and r analogous to the usual volume deceleration parameter and 
“Hubble spatial gradient” as follows 


d 0 /3 := ~{l + q)P, di/3:=-r/3. 


(3.10) 


so that 

doEi 1 = E\ l (l + q) + ^ 2 -eo(ei 1 ), doQ = <2(1 + q) + ~^~p e o(&)> 

d 0 E = S(1 + q) + -jpeo(a+), 8 qU = U(1 + q) + —^e o(h), 8qA = A{ 1 + q) + -^eo(a), 

dolC = 2/C(l + q) + ^g 3 e o{K)i d 0 n = 217(1 + q) + ^o e o (A), ®o v = 


and 

0iQ = 'rQ + ^2 e i( 0 )’ = rS + (^2 e i( cr +)’ 01^ = rU + -^e 1 (u), diA = rA + -^ei(a), 

0!/C = 2r/C + ^ 3 -ei(X), 0il7 = 2rl7 + 0,'u = iv. 
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Thus, the field equations reduce to 


d 0 El = (q + 3Z)El, 
d 0 IC = 2qK, 

■ (3cg + 2C(j) d\Q ■ r(3cg + 2c a )Q 3 7 c a v£t 

c a d 0 U -7777-77- = c a U(q - 1)-77;-77-b 


(2c* - 1) 


(2c ct -1) (2c ct — 1) (1 — v 2 ) ’ 


(3.13a) 

(3.13b) 

(3.13c) 


d 0 Q _ ^+; X ) = g( i + q _ Q) _ ^ + W _ 2 !) ^ + 


3 (3c e + 1) 


3 (3c 0 + 1) 3 (3cg + 1) 


| (c a + 1 ) U 2 | 2 (2c ct — 1) S 2 | 12(( 7 -2)u 2 -3 7 + 2) 


dn£- 


3 (3 Cg + 1) 3 Cg + 1 

(c a + l)d!& 


2 (3c e + 1) (1 — V 2 ) 


2 \ ’ 


3(2 C(T - 1) 


= £(! + ,- 3g) - r ±“ + V? + ( 1 „I 2C “ M . C/ + 


+ 


3/C 


+ 


(5c a + 2) t / 2 


3(2c ff - 1) 3(2co- - 1) 

-4 2 | (3cg + 1) Q 2 | 0((7 + 1)u 2 -1) 


2 ( 2 co- — 1 ) 6 ( 2 ^ — 1 ) 2 ( 2 co- — 1 ) 2 ( 2 c* - 1 ) 2 ( 2 c* - 1 ) (1 - u 2 ) 


d 0 A + (3<=7+2^)9 iC = qA + r (3c, + 2c„) Q 


37^12 


2c* - 1 


2Ccr — 1 2 (2c* - 1) (1 - V 2 ) 


-U, 


- (7 — 2)uc/il2 'yQdiv 2 r yAvCl 2qCl (l — (7 — l)v 2 ) 

c/o$i H-77-77- =-77-1-77-b 


G- 


G _ 




G _ 


3 7 Q (l — u 2 ) 12 2(7 — 2)rul2 ^ 212 ((1 — 2 7 )u 2 + l) 


G_ 




G _ 


5 0 u - 


(7 - 1 ) (l - u 2 ) 2 5 i !2 (7 - 2 )u 5 iu _ 2(7 - 1 )A (l - v 2 ) 


7 12G_ 


+ 


G_ 


3(7 — 2 )Q (l — v 2 ) v 2(7 — l)r (l — v 2 ) 2 (l — u 2 ) v . ( 

- gZ - W- - — 72 - + U{l-v ) 


G- 
v 

gZ 


(3.13d) 


+ -£ 2 
2 ’ 

(3.13e) 

(3.13f) 


(3.13g) 


(3.13h) 


and 


dxN~ l = (r-U) AC" 1 , 
dxlC = 2(r + A)JC, 

d\A + c a d\U = — —JC + r (c a U + + 2c a AU — —c a U 2 + —A~-\- 

3 , ~ ^ 312 ((7 - l)u 2 + 1) 3 

— o (3c$ + 2c*) Q" + 3 (2c* — 1) Q H- - -77-3c* + —, 

z z (1 — v z ) z 

8l £ - ( 3 7+ 2) 3l g = rS + MS - r (3 7 + V Q + 3lVSl 


2 c* — 1 
doA = (q + 3£)A. — U + r. 


2(7 - 1 


2 ( 2 c* — 1 ) (1 — u 2 ) ’ 


(3.14a) 

(3.14b) 


(3.14c) 

(3.14d) 

(3.14e) 


Note that A 2 ,1C, only appear in the equations (3.13d), (3.13e), (3.13f) and (3.14c) via 
the combination T> = .A 2 — 3/C (and in the other equations also via terms of the form AU, Av, 
the definition of q (B.l), etc.). We further develop the governing equations; however, since 
this is rather technical we continue this development in Appendix B. 
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4 Dust models 


Let us consider dust models with 7 = 1 (p = 0). The special case of dust, in which the 
governing equations simplify considerably, is of particular interest at late times. In GR we 
immediately obtain the simple Lemaitre-Tolman-Bondi (LTB) model with u = 0 (and since 
N = N(t) is a function of t, we can set N = 1 by a time rescaling), where the fluid is 
“comoving” (v = 0); i.e., u = 0 and v = 0 simultaneously (see Appendix C). This is not 
possible in the models here; in general v cannot be zero in the dust case and there is no 
GR-like “LTB” model. 

If v = 0 (7 = 1), from equations (3.5c) -(3.5h) we immediately find that u = 0, which is 
a contradiction (for equations (3.5c) -(3.5h)). Therefore, in our formalism, the perfect fluid 
must be tilting (u / 0). In general, we thus need to investigate dust with u/0 and u A 0 
(i.e., non-“LTB”). Let us next consider the case u = 0 (see equations (B.21a)-(B.21g)). From 
equations (B.22a)-(B.22d), if v = 0 we then find that (either) 3 cq + 2 c a = 0 (eg = c a = 0 in 
GR) (or ei(0) = 0, which is valid in the spatially homogeneous models); we could investigate 
this special model further. 

Let us also consider the subcase U = 0 and v = 0 in normalized variables. 

4.1 Normalized equations 

Let us study the special subset U = v = 0, with 7 = 1 (see Appendix B). We also assume 
that cq = 0 and c a / 0 : 


d 0 IC = 2 qK, 

doQ = Q(1 + q ~ Q) + 2 {2c a - 1) (1 - Qf - 


(4.1a) 

(4.1b) 

(4.1c) 
(4-Id) 


do A = qA, 

d 0 n = (2q-3Q + 2)n, 


subject to the restrictions: 


diK = 2{r + A)IC, 

d\A = ~~7}ie + i"A + ~^A~ — 3C(j<2" + — (2co- — 1 ) (2 Q — 1 ) + —12, 
d\Q = rQ , 


(4.2a) 

(4.2b) 

(4.2c) 


where q and r are defined by: 



(4.3) 


r = -3A(l-Q). (4.4) 


(4.4) 


We also have that do ■= A f 1 dt, d\ := E\d x . [The only remaining freedom is the coordinate 
rescalings t —>• f(t) and x —>• g(x)\. 
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4.2 Special dust model 

The terms 1C, A 2 , only appear in the evolution equations for doSl,doQ via (through q) the 
combination V = A 2 — 3/C. Hence (assuming c CT A 0) we have the reduced (closed) evolution 
system: 

d 0 V = 2 qV, (4.5a) 

d 0 Q = Q(1 + q - Q) + 2 (2- 1) (1 - Q) 2 - h, (4.5b) 

d 0 Cl = (2 q -3 Q + 2)Sl, (4.5c) 

where 

q = 1 A v-2c a ( 8 c ff - 3) Q 2 + \Qc a (2ca - 1 )Q 

2 ( 2 c a - 1 ) f 

+ 2c a O + (1 — 8c a )(2c a — 1)|. (4.6) 

In the decoupled evolution equations above (which are only valid strictly speaking for 
c a A 0) we have not yet applied any constraints. The constraint eqns for “LTB”-like models 
(i.e., dust models in Einstein-aether theory with U = 0 and v = 0) imply either 3cg + 2c a = 0 
or ei(0) = 0. The problems regarding “LTB” come from the constraints for ei(0) and ei(<7+); 
when normalizing with 3/3 = @ + 3cr+, these constraints get hidden in the normalized variables 
(because 6 decouples in the normalized eqns and £ is related to cr + /{0 + 3<r+)), and so there 
are no problems per se with normalized equations But they do not represent any “LTB” 
model because the constraints are not satisfied. 

The FLRW models in this special dust model (with c a A 0) have V = 0,Q = 1, Q = 1 
(q = \) and correspond to an equilibrium point (the point P 4 below). In the Kantowski- 
Sachs models, A = 0, r = 0 , and there are no spatial derivatives, and the constraints can 
be used to eliminate the (non-zero) T> and the resulting system becomes 2 -dimensional (the 
Kantowski-Sachs models will be studied later using a different normalization). 

Assuming A A 0 , we can define the new spatial derivative d v = A~ 1 d\, whence the 
spatial derivatives become: 

d n V = 3 (— 2c ct (1 - Q ) 2 + {V - 1)(-1 + 2 Q) + Si) , (4.7a) 

d v Q = —3(1 — Q)Q. (4.7b) 

The commutator equation is given by 

[d v ,d T \ = A _1 (qd 0 + [9i,a 0 ]). (4.8) 

There is no spatial restriction for SI] thus, it is freely specified at the initial spatial 
hypersurface. 

Let us summarize the equilibrium points of the system (4.5) and their eigenvalues (see 
Table 1), and discuss their stability: 

1. Point Pi exists (i.e., with — 1 < Q < 1) for c a < 0 or 0 < c a < A It is a source 
for c a < 0 or 0 < c a < |; a saddle for | < c a < (The equilibrium points are 
non-hyperbolic for other values of the parameter c a )- 
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Table 1. Equilibrium points and eigenvalues of the system ( 4 . 5 ). 


2. Point P 2 exists for | < c a < It is a sink for | < c a < a saddle for ^ < c a < ^. 

3. Point P 3 exists for c a > It is a sink for c c > a saddle for < c a < 

4. The [FLRW] point P 4 always exists and it is a saddle. 

5. The point P 5 always exists and it is a saddle for c a / 0. 

6 . The point Pq always exists. It is sink for c a <\ [two complex conjugate eigenvalues 

with negative real part for 0 < c a < ^]. It is a saddle for c a > \. For c a = \ it is a 
saddle too. 

7. The point P^ exists for c a < |. It is a source for | < c a < |. Non-hyperbolic for 
c a € {0, |}. Saddle otherwise. 

Discussion: Let us define 

D a = ( 1 - 8 c cr )(l - 2 Ca) + 16c*(1 - 2co-)Q - 2c a (3 - 8 c a ) Q 2 , 
and consider the equilibrium points at finite values: 

(&) D = 0,D = 0, Q = Q a (where there are constraints on the parameter c a in order for Q a 

to be physical) [the points P\,P 2 ,P 3 in Table 1]; generically a saddle. 

(b) V = 0, n = 1, Q = 1 (q = ^). Eigenvalues: Ai = 1, A 2 = 1 and A 3 = — | The FLRW 
equilibium point [point P 4 in Table 1] is always a saddle. 

(c) V = 1 — 2co-, fl = | c a , Q = | (q = 0) [point P 5 in Table 1.] 

(d) V = D a ,Q = 0 (q = |), and either (i) Q = 1 (T> = 1; no shear) or (ii) Q = Qd = 
2(1 — 2c CT )/(3 — 4c,j) (7? = 3(1 — 2c ff )(3 — 8c fr )/(3 — 4c,j) 2 ). Eigenvalues: (i) (Q = 1 [point Pq 

in Table 1]) Ai = —1, A 2,3 = —1 ± J j (negative real part for c a < ^) - corresponding 
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Figure 1. Phase space of the system ( 4 . 5 ) for the choice c a = The sinks are P2 and Pq. P7 and 
Pi coincide; they are non-hyperbolic and behave as the sources. 



Figure 2. Phase space of the system ( 4 . 5 ) for the choice c 0 = 0.45 £ (|, )r§). The sinks are P2 and 
Pq, and P7 is the source. 


to a sink, (ii) [point P 7 in Table 1] Ai = y, A2 = —1 + 0 (c a ), A3 = 1 + 0{c a ) - which 

is a saddle for small c a . 

Summary of sinks: Pq for c a < P2 for | < c a < ^|, P3 for c a > ^|. In all cases 
—>• 0 to the future. For P 2 and P 3 , V —>• 0, but for Pq, —>• 1 CQ —>• 1^ and the shear goes 
to zero at late times (for small c a ). There is a range of values of the parameter c a for which 
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Figure 3. Phase space of the system (4.5) for the choice c a = 0.7. The sink is P 3 . 


the sinks P2 and P3 represent inflationary solutions. 

For illustration, we present some phase portraits of the system (4.5) for the x-constant 
surfaces in figures 1-3. In figure 1 we present some orbits of the phase space for the parameter 
c a = |. The sinks are P2 and Pq. P7 and P\ coincide; they are non-hyperbolic and behave as 
the sources. In figure 2 we present the evolution of the system (4.5) for c a = 0.45 € (|, g|). 
The sinks are P2 and Pq. T 7 is the source. Finally, in figure 3 we present the phase portrait 
for the choice c a = 0.7. The sink is P3. 

5 Special cases with extra Killing vectors 

Spherically symmetric models with more than 3 KVF are either spatially homogeneous or 
static. Spatially homogeneous spherically symmetric models are either Kantowski-Sachs mod¬ 
els, or the Friedmann-Lemaitre-Robertson-Walker (FLRW) models (with or without cosmo¬ 
logical constant A) [or locally rotationally symmetric (LRS) Bianchi I and Bianchi III models]. 
For the FLRW and Kantowski-Sachs models we use the equations in the case that u = 0, 
which follows immediately from the condition that N = N(t). Static and self-similar spheri¬ 
cally symmetric models have been studied in [21, 46]. 6 

5.1 The FLRW models 

For the FLRW models the source must be of the form of a comoving perfect fluid (or vacuum) 
and the aether must be comoving. The metric has the form 

ds 2 = — N(t) 2 dt 2 + l 2 (t)dx 2 + l 2 (t) f 2 (x){d‘d 2 + sin 2 i!)d(p 2 ), (5.1) 

with 

f(x) = sinx, x, sinhx, (5-2) 

6 Recall that we have chosen a gauge so that the aether is aligned with eo. 
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for closed, flat, and open FLRW models, respectively. The frame coefficients are given by 
ei 1 = and e 2 2 = £~ 1 (t)f~ 1 (x). Then a + = ^eo ln(ei 1 /e 2 2 ) vanishes. Furthermore, 

a = — ■ N = N(t) implies that u = 0; i.e., the temporal gauge is synchronous, and 

we can set N to any positive function of t (we usually choose N = 1). The Hubble scalar 
H = e o \n£(t) is also a function of t. ‘ 

For the spatial curvatures, 3 S+ does vanish because (5.2) implies eia = AT , 8 while 3 i? 
simplifies to 

3 R = _, « = 1,0,-1, (5.3) 

for closed, flat, and open FLRW, respectively. The evolution equation for er + and the Codazzi 
constraint then imply that 7 T+ = 0 = q\. 

FLRW cosmological models with aether and a comoving perfect fluid comoving have 
been studied previously [1, 3, 6 , 37]. It was found that there is no essential affects on 
standard cosmology in the minimal aether theory. FLRW cosmological models with a scalar 
field were studied in [2, 16, 17]. The decay of tilt has also been studied in (anisotropic and 
non-comoving) models with A [4, 5]. 

5.1.1 The FLRW models in normalized coordinates 

The FLRW models in normalized coordinates are characterized by £ = 0 (Q = 1 ), U = 0, 
v = 0, and we can use the remaining coordinate freedom to set N = 1 (where f3 = -^). 
We recall that d\ := E\d x and doE\ = qE\. We use the remaining spatial freedom to 
simplify f(x) as in equation (5.2) for the FLRW metric as above, where (d x f(x )) 2 = 1 — k/ 2 , 
d x d x f{x) = —nf, and so we obtain: 9 


1C 


A 


N 2 

3 i 2 f 2 

(5.4) 

Nd x f 
if ’ 

(5.5) 


and hence 


A 2 - 3JC = diA = ^-K. (5.6) 

£2 ’ N 

5.1.2 The subset U = v = 0 with Q = 1 

We take the equations in the case U = v = 0 presented earlier, and set Q = 1. We again 
assume that eg = 0 and c a ^ 0 (and, in principle, 7 / 1). Since Q = 1, £ = 0 (i.e., the shear 
is zero), which is not in general an invariant set. We immediately have that r = 0, whence 
diM ^ 1 = 0, and we can rescale time so that N = 1 and do := d T , where r is essentially 
logarithmic time. We also have that Cl is independent of space, and that 

7 We shall not list the KVFs as they are complicated in spherically symmetric coordinates and not needed 
here. 

s That ei a does not vanish is consistent with the frame vector ei not being group-invariant. 

9 Flat FLRW power law models: In the flat case k = 0, f(x) = x (and r = InR)). At the equilibrium points 
we have that l = t p . 
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subject to the restrictions: 


d T El = qE {, 

(5.7a) 

d T IC = 2 qK, 

(5.7b) 

d r A = 

(5.7c) 

3 r 12 = (2g - 37 + 2)12, 

(5.7d) 


diK = 2AJC, (5.8a) 

dl A = -^IC + ^A 2 + ^(n-l), (5.8b) 

(5.8c) 

where q is defined by: 

q = h (3 7 - 2), (5.9) 

and where 

0 = -3JC + A 2 + n- 1. (5.10) 

Note that if we differentiate this constraint and use the constraint and the definition 
of q, we obtain zero; hence the constraint is conserved along the evolution. We can use 
this constraint to eliminate /C from the above equations. We note, as expected, that all 
dependence on c a has dropped out. We also note that the above system has the equilibrium 
points n = 0, q = 0, corresponding to late time vacuum, and Q = 1, q = ^(3q — 2), the early 
time flat solution. 

Finally, using equations (5.4, 5.5 5.6), we obtain 

d r Q = (2q - 3q + 2)fi, (5.11) 

where q = (37 — 2) and Cl = 1 + 2|, as expected. 

Note that in the general solution (i.e., not FLRW) we can define T> = A 2 — SIC, so that 
d T V = 2 qT>, subject to the restrictions: d n T> = 3 (T> + Cl — 1), where we have introduced the 
new spatial coordinate d n = d\/ A. Since 0 = 11 + 12 — 1 we then obtain d ri T> = 0, which 
implies d n il = 0. We then obtain 

e 2r _ (37 - 2)e 2r 

e 37T+ci _|_ g 2 t ’ ^ 2 ( e 37' r +Ci g2r) ’ 

e 37r ( 02 ( 77 ) 2 — e Cl ) 

3 (g 37 T + c i + e 2r ) 

The equations 

3,,/C = 2/C, (5.12a) 

Ad v A = h>C + ^A 2 + ^{Cl-l), (5.12b) 


V = 


1 


3—37T—ci+2r _j_ ^ 


n = 


«- , C2{V) ,K = 

V e 3 7 i"+ci _|_ g 2 r 
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are identically satisfied if 


e Cl + c 2 (r]) (c' 2 ( 77 ) - 02 ( 77 )) = 0 , c 2 (rj) / 0. 


(5.13) 


The above equations admit the solutions 

oM = e; = 

p37r+2c2+2r; \/p c l p2(c2+?7)p^2“ 

jr = _4 = 4 - v - _ 

3 ( e 3 7 r+ci + e 2r) ’ y^r+ci + g 2r 


5.2 The Kantowski-Sachs models 


We now investigate the spatially homogeneous subcase, in which a full global analysis is pos¬ 
sible. It is of particular interest whether general solutions can asymptote towards spatially 
homogeneous solutions at late or early times. The spatially homogeneous spherically symmet¬ 
ric models (that has 4 Killing vectors, the fourth being d x ) are the so-called Kantowski-Sachs 
models [36]. We shall consider the special comoving aether case. The metric (2.1) simplifies 
to 

ds 2 = — N(t) 2 dt 2 + (ei 1 ^)) -2 ^ 2 + (e 2 2 (t ))~ 2 (dd 2 + sin 2 dd<p 2 ); (5-14) 

i.e., N, ei 1 and e 2 2 are now independent of x. The spatial derivative terms ei( ) vanish and 
as a result a = 0 = u. Since u = 0, N is a positive function of t which under a time rescaling 
can be set to one. This metric choice forces the fluid to be non-tilted (v = 0) [assuming 
(*> 0,7 > 0]. 

The evolution equations for the Kantowski-Sachs metric for an Einstein-aether spheri¬ 
cally symmetric cosmology, in the presence of a perfect fluid, are: 


eo(ei 1 ) 
eo(K) = 

eo( 0 ) = 

eo(cr+) 
e o(A) = 


= -§(0 - 6u + )ei 1 
= -§(0 + 3 a+)K 
e 2 6(2 C(T - 1)cj 2 

~~3 + 3 c e + 1 


+ 


(2 ~ 3 7 ) A 
2(3c0 + 1) ’ 


= (3 c e + 1)0 2 

9(2c ff - 1) 
-7 


- 6 o+ - cr 2 + 


A 4 


3(2 C(T - 


1 )’ 


(5.15a) 

(5.15b) 

(5.15c) 

(5.15d) 

(5.15e) 


with the constraint 


k+ (3c + 1 = a _ 3(2Cj _ !)4. 


(5.16) 


We choose the following normalized variables (which are bounded for 1 — 2co- > 0; note 
that we do not use the /?—normalization for convenience here): 




x = ^y = 


V3a + _ y/K ^ _ 9 

I) ,Z D M ^3 D 


(5.17) 


where 


D = 



(5.18) 
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and the new time variable f = -^eo(/). 

We then obtain the full 4 dimensional (4D) system: 

x , = x ( Qy 2 (2 - 4 c a ) + Q^_ + Q ( 2z2 ~ 3 t) y-Y (37 - 2)Qx 3 

X ^ \/3 (3c 0 + 1) \/3 2\/3 V3/ 2x/3(3c e + l) 

, _ V3Q 2 (—3cg - 2y 2 c ff + y 2 - 1) 2 / (37 - 2)Qg/ _ 1 

V ~ 6co- — 3 X V2\/3(3c e + l) + V3(2c a -1) 

2\[3Qy (3 c e + 2 y 2 c a - y 2 + l) 

9c0 + 3 

^ _ Y ^VSQy 2 (1 - 2cp-) _ Q^y\ + (37 - 2)Qx 2 z 
^ _ ~ 9c e + 3 V3; + 2\/3 (3c 0 + 1) 

n , = 2 ( 2\/3y 2 (2cp- - 1) Qy\ (2 - 3j)x 2 z 2 
y 9c 0 + 3 >/3/ 2\/3 (3c 0 + 1) ’ 

The variables (5.17) are related through the constraints 

-3c e Q 2 + x 2 - (2 c a - 1 )y 2 = 1, 

Q 2 + Z 2 = 1, 


(5.19a) 

(5.19b) 

(5.19c) 

(5.19d) 

(5.20a) 

(5.20b) 


which are preserved by the 4D system. From the equations (5.20) it follows that Q and 2 are 
bounded in the intervals Q e [—1,1], 2 € [0,1] (for expanding universes Q > 0). However, 
since 1 — 2c CT is not necessarily non-negative it follows that x and y are unbounded, unless 
1 - 2 c a > 0. 

The restrictions (5.20) allow the elimination of two variables, say x and z. This leads 
to the following 2-dimensional dynamical system: 


, _ V3Q 2 V3Qy (eg (( 3 7 — 2)Q 2 — 4) + 7 — 2) ^3(7 - 2)Qy 3 (2c a 

^ 3 — 6 c a 6 eg + 2 6 cq + 2 

1 _ (■ Q 2 - 1) y 2 

+ x/3 (2c a - 1) y/3 

= V3(2 - 37) ^(37 - 2 )Q 4 c d _ V3(3 7 - 2)Q 2 (3 c e - 1 ) 

4 18 c 0 + 6 6c e + 2 18 c e + 6 

V^( 7 ~~ 2 ) (Q 2 — l) y 2 (2c a — 1) (1 -Q 2 )Qy 

6 Cq + 2 y /3 


(5.21a) 


(5.21b) 


We shall study the general case in future work (using the /3— normalization). Let us 
consider the following special case here. 


5.2.1 Special case. 

Let us assume 

3 c e = ci + 3c 2 + c 3 = 0 (5.22) 

(see Appendix A and the references [16-18]), and define c 2 = 1 — 2 c„ = 1 — 2(c\ + c 3 ) > 0. 
This choice leads to a compact phase space. 
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With these special values of the c’s, the evolution equations for Kantowski-Sachs models 
simplify and the constraint becomes 


9 2 2 2 

K -\—— — fi + 3c 
o 


The following normalized variable 


V 1 = 


\/3ccj + 


L> 


is chosen for convenience, whence the variables are related through the constraints 

x 2 + y\ = 1, 

Q 2 + z 2 = 1. 


(5.23) 

(5.24) 

(5.25a) 

(5.25b) 


Thus, the phase space is compact with x £ [—1,1], y\ £ [—1,1] and Q € [—1,1], z 6 [0,1] 
(for expanding universes Q > 0). 

The system for (yi,Q) reduces to 


Vi = 


Q' = - 


{yj - 1) (3c(7 - 2)Qyi + 2Q 2 - 2) 

2y/3c 

(' Q 2 ~ 1) (c (-37 + 3(7 ~ %i + 2) + 2Qj/i) 
2 V^c 


(5.26a) 

(5.26b) 


Since the evolution equations are invariant under the transformation y\ —>• — yi and 
c —>• —c, without loss of generality we can assume c > 0. Scaling the time derivative by the 
positive factor 2\/3c, we then obtain: 


2/1 = - (2/? - 1) (3c( 7 - 2)Q yi + 2 Q 2 - 2), (5.27a) 

Q' = ~ (Q 2 - 1) (c (-37 + 3(7 - 2 )y 2 + 2) + 2 Q Vl ) (5.27b) 


In tables 2 and 3 we present the equilibrium points of the system (5.27) and discuss 
their stability. We have that c> 0,-1 < Q < 1,-1 < y± < 1. Some of the equilibrium 
points do not exist for certain values of ”c”. We have not analyzed the non-hyperbolic ’’stiff 
fluid” case, 7 = 2, in which there are zero eigenvalues. Clearly, the case c = 0 is not included 
here (the GR case), since the equations are not valid in that case. 

Let us enumerate the stability conditions for the hyperbolic equilibrium points: 


1. The equilibrium point Pi is a source for c > 0,0 < 7 < |, and a saddle for | < 7 < 2. 
Non-hyperbolic for 7 = | or 7 = 2. 

2. The equilibrium point P 2 is a sink for c > 0, 0 < 7 < |, and a saddle for | < 7 < 2. 
Non-hyperbolic for 7 = | or 7 = 2. 

3. The equilibrium point P3 is a sink for 0<7<2, c > ^, and non-hyperbolic for c = ^ 
or 7 = 2. Saddle otherwise. 

4. The equilibrium point P 4 is a source for c > 0,0 < 7 < 2. Non-hyperbolic for 7 = 2. 

5. The equilibrium point P5 is a sink for c > 0,0 < 7 < 2. Non-hyperbolic for 7 = 2. 
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Label 

Coordinates: (yi,Q) 

Eigenvalues 

Pi 

( 0 ,- 1 ) 

3c(2 — 7 ), 2c(2 — 3 7 ) 

P 2 

( 0 , 1 ) 

1 

CO 

0 

1 o 

1 

1 

to 

C> 

"to 

1 

CO 

P 3 

(- 1 ,- 1 ) 

—6c(2 — 7), 4(1 — 2c) 

Pi 

(- 1 , 1 ) 

6c(2 — 7), 4(1 + 2c) 

P 5 

( 1 ,- 1 ) 

—6c(2 — 7), —4(1 + 2c) 

P 6 

( 1 , 1 ) 

6c(2 — 7 ), —4(1 - 2c) 

Pr 

(— 1 , — 2 c) 

2(4c 2 - 1),4 [c 2 ( 37 - 2) - l] 

Ps 

( 1 , 2 c) 

—2(4c 2 - 1), -4 [c 2 (37 - 2) - l] 

P 9 

M 2 - 37 ) 2 A 

y d ’ dj 

c(—e— 37 + 6 ) c(e—3 7 +6) 
d ’ d 

P 10 

( c ( 2 - 3 7 ) 2 \ 

V d ' d) 

c(—e+ 37 — 6 ) c(e+3 7 — 6 ) 
d ’ d 


Table 2. Equilibrium points of the system (5.27) and their eigenvalues. We use the notation d = 
^3(7- 2)( 37 - 2)c 2 +4 and e = \/{2 - 3 7) 2 - 27 7 + 22 - 

6 . The equilibrium point P$ is a source for 0 < 7 < 2 ,c> 4 .A saddle for 0 < 7 < 2,0 < 
c < -j. Non-hyperbolic for 7 = 2 or c = |. 

7. The equilibrium point P7 exist for 0 < 7 < 2 , 0 <c<^. It is a sink for 0 < 7 < 2, 0 < 
c < Non-hyperbolic otherwise. 

8 . The equilibrium point P$ exists for 0 < 7 < 2 ,0 < c < 1. It is a source for 0 < 7 < 
2, 0 < c < Non-hyperbolic otherwise. 

9. The equilibrium point P 9 exists for 0 < c < ^,0 < 7 < |, or 0 < c < 5,7 = 2, or 

c>^, 0 < 7 <|. It is a saddle for 0 < 7 < |,c > 0. Non-hyperbolic for 7 = | or 

7 = 2 . 

10. The equilibrium point P 10 exists for 0 < c < ^,0 < 7 < or 0 < c < ^,7 = 2, or 

c > ^,0 < 7 < It is a saddle for 0 < 7 < |,c > 0. Non-hyperbolic for 7 = | or 

7 = 2 / 

Discussion. In the case c a < ^ (i.e., c > 0 ), when 7 < P 2 is the unique shear-free, 
zero curvature (FLRW) inflationary future attractor, and for | < c a < \ (be., 0 < c < |) 
and 0 < 7 < 2 the sources and sinks are, respectively, P4 & Ps and P5 & P 7 . All of these 
sources and sinks have maximal shearing and all, except P7, have zero curvature; the sink P7 
does not have zero curvature. For c a < | (i.e., c > |) the points P7 & P$ do not exist, and 
the sources and sinks with maximal shearing are P4 and P5, respectively. 

In figures 4 - 7 we present some orbits in the phase plane of the system (5.27) for 
different choices of the parameters. In figure 4, 7 = 0 and c a = 0.3. The sinks are P 2 ,P 5 
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Eq. Pt. 

7 value 

0 < c < \ 

c=\ 

\<c 

Pi 

0 — 7 < f 

++ source 

++ source 

++ source 


II 

COllsJ 

+0 

+0 

+0 


| < 7 < 2 

H— saddle 

H— saddle 

H— saddle 


7 = 2 

0 - 

0 - 

0 - 

P 2 

0 — 7 < | 

- sink 

- sink 

- sink 


-3 

II 

CO| LO 

-0 

-0 

-0 


| <7<2 

— |- saddle 

— b saddle 

— b saddle 


7 = 2 

0 + 

0 + 

0 + 

P3 

0 < 7 < 2 

—|- saddle 

-0 

-sink 


7 = 2 

0 + 

00 

0 - 

P 4 

0 < 7 < 2 

++ source 

++ source 

++ source 


7 = 2 

0 + 

0 + 

0 + 

P5 

0 < 7 < 2 

- sink 

- sink 

- sink 


7 = 2 

0 - 

0 - 

0 + 

P 6 

0 < 7 < 2 

— bsaddle 

+0 

++ source 


7 = 2 

0 - 

00 

0 + 

Pi 

0 < 7 < 2 

-sink 

-0 

DNE 


7 = 2 

-sink 

00 

DNE 

P 8 

0 < 7 < 2 

++ source 

+0 

DNE 


7 = 2 

++ source 

00 

DNE 

P 9 

0 — 7 < f 

—b saddle 

—bsaddle 

—bsaddle 


-3 

II 

CO|K 3 

+0 

+0 

+0 


7 = 2 

00 

00 

00 

Pi 0 

0 — 7 < | 

—b saddle 

—b saddle 

—b saddle 


2 

1 3 

-0 

-0 

-0 


7 = 2 

00 

00 

00 


Table 3. Stability of the equilibrium points of the system (5.27). 


and P7. The sources are P\,P± and Pg. P 3 ,P&,Pg and P10 are saddles. In figure 5, 7 = 0 
and c a = 0.6. The sinks are Pg. Pg and P5. The sources are P\ , P4 and Pq. P7 and Ps do 
not exist. The saddles are Pg and Pro- In figure 6 we present the phase plane of the system 
(5.27) for the choice of parameters 7=1 and c a = 0.2. The sinks are P5 and P7. The sources 
are P4 and Pg. The saddles are Pi, P2, P3 and Pq. Pg and P10 do not exist. Finally, in figure 
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Figure 4. Phase plane of the system ( 5 . 27 ) for the choice of parameters 7 = 0 and c a = 0 . 3 . The 
sinks are Pj, P5 and P7. The sources are Pi, Pa and P$. P3,P e ,P 9 and Pio are saddles. 



Figure 5. Phase plane of the system ( 5 . 27 ) for the choice of parameters 7 = 0 and c a = 0 . 6 . The 
sinks are P2,P?, and P5. The sources are Pi, Pa and Pq. P7 and P§ do not exist. The saddles are P9 
and Pio. 

7 , 7 = 1 and c a = 0 . 6 . The sinks are P3 and P5. The sources are P4 and Pq. Pi and P2 are 
saddles. The points P7-P10 do not exist. 
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Figure 6. Phase plane of the system (5.27) for the choice of parameters 7 = 1 and c a = 0.2. The 
sinks are P5 and P7. The sources are P4 and Pg. The saddles are Pi,P 2 ,P3 and P§. Pg and P10 do 
not exist. 



Figure 7. Phase plane of the system (5.27) for the choice of parameters 7 = 1 and c a = 0.6. The 
sinks are P3 and P5. The sources are P4 and Pq. Pi and P 2 are saddles. The points P7-P10 do not 
exist. 


6 Static models 

Models that are not evolving with time are also of physical importance, although perhaps 
more from the astrophysical point of view than from the cosmological one. In particular, 
much physical information can be obtained from a qualitative analysis of the models. Let us 
consider the static case eo(-) = 0, for a mixture of a perfect fluid and a scalar field. [In this case 
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u = din N/dx ^ 0 and the perfect fluid is forced to be non-tilted (v = 0).] We will consider 
barotropic equations of state jj, = Since 9 = <j+ = 0 in the static subcase, V depends 

only on the scalar field <f> = <f>(x). Furthermore, the irreducible components of the scalar field 
energy-momentum tensor are given by p? = ^ei(^) 2 + V,p^ = — |ei (<)>) 2 — V,q^ = 0 and 
*+ = -|ei(0) 2 [15]. 

The equations for the variables a, u,p, 4>, K, N are: 

ei ( a ) = of + 7TT + e i(^) 2 -r + 2c a u 2 + 3 ail + K, (6.1a) 

2(c a + l) c a + l 

e i («) = of + -+ 2au-u 2 , (6.1b) 

2 (c a + 1) c a + 1 

e i (£) = -<i(A + £) (6.1c) 

ei(ei(0)) = - (it - 2a) ei(</>) + V^, (6.Id) 

ei(/\) = 2aX, (6.1e) 


where V^, denotes differentiation with respect to cj). The system satisfies the restriction 

a 2 = c a ii 2 + 2ait + p + ^ei(^) 2 — V + K. (6.2) 

Taking the differential operator ei(...) of both sides of (6.2), using the equations (6.1) to 
substitute for the spatial derivatives, and again using the restriction (6.2) solved for K , we 
obtain an identity. Thus, the Gauss constraint is a first integral of the system. The aether 
constraint is identically zero. 

Let us now show how equations (6.1) can be used to obtain exact solutions and, ad¬ 
ditionally, use the dynamical systems approach to investigate the structure of the whole 
solution space. The first thing to do is to select a suitable radial coordinate. We may choose 
a new radial coordinate A such that the equation (6.1e) has a trivial solution. A reasonable 
way to do this is to select an r-coordinate such that ei(/) = — ard r (f) (as in [47]). This 
implies K oc r -2 and f ei 1 (x)~ 1 dx = — f a _1 dlnr. As we will see later, for the dynamical 
systems investigation it is better to use the new time variable r = lnr, which takes values 
over the whole real line. 

Here we shall study the two special cases: (i) perfect fluid (previous work has assumed 
a comoving aether and a comoving fluid), (ii) vacuum (stationary with a scalar field and 
harmonic potential). In particular, for the stationary aether case, it is also of interest to 
choose a frame in which the aether is non-comoving. For a non-comoving stationary aether 
and a tilted fluid, it follows [in an analogous way to the static case] that the perfect fluid 
must be non-tilted (v = 0). Additionally, since 6 = = 0 in the stationary subcase, V 

depends only on the scalar field <f> = (f>(x). We shall present a more comprehensive analysis in 
[46]; in particular, the “evolution” equations for the tilt a, a and ii are given in the Appendix 
therein. 


6.1 Static case with perfect fluid with linear equation of state and no scalar 
field. 

To investigate this model we will use the approach of [48]. First, let us consider no scalar 
field in (6.1) and use the linear equation of state 

A = Mo + (r? - 1 )P, (6.3) 
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where the constants hq and 77 satisfy /oo > 0 ,77 > 1. The case 77 = 1 corresponds to an 
incompressible fluid with constant energy density, while the case //o = 0 describes a scale- 
invariant equation of state. 

Introducing the new dimensionless variables 


Mo u 

X\ 9 , %2 5 

a 


X 3 = 


P_ 

,2’ 


X4 = 


K 

l 2 > 


we obtain the dynamical system 
dx\ 


dr 


= x 1 


Xl +x 3 (ji + 2) + 4x 2 Ca + 2 X4 + 6x2 ) , 
C a + 1 


dx 2 _ (x 2 - l)(xi + 0:3(77 + 2)) , 0 3 


dr 2 (c a + 1 ) 

dx 3 _ x 3 (xi + 0:3(77 + 2)) 
dr c a + 1 

dx 4 


T 2 x 2 C a + 0:2(074 + 4 x 2 — 2 ), 


+ 4x 3 X2C a + x 2 (xi + x 3 (?7 + 6)) + 2x4X3, 


dr 


= X 4 


x\ + x 3 (?7 + 2 3 4 ) + 4^2^ + 2x 4 + 6x 2 - 2 1 , 


C a + 1 


subject to the constraint 


1 = X2Q1 + 2x 2 + 0:4 + x 3 . 


(6.4) 


(6.5a) 

(6.5b) 

(6.5c) 

(6.5d) 


( 6 . 6 ) 


The constraint ( 6 . 6 ) is preserved by the dynamical system (6.5). Solving the constraint ( 6 . 6 ) 
for x 4 and substituting back into the system (6.5) we obtain the reduced system: 


dx 1 X1X3 (77 — 2 c a ) 9 

1 = 3VI a) + 2xiX 2 Ca + Xl 

dr c a + 1 


Xl 


C a + 1 


T2 ) +2 xix 2 , 


dx 2 
dr 
dx3 
dr 


= o : 2 


= 0:3 


Xl 


2 c a + 2 

Xl 


+ 0:3 


77 + 2 


2 c a + 2 


-1 -1 - 


07l X 3 (77 + 2) 


2 c a + 2 2 c a + 2 


(6.7a) 


T xlc a + 2x|, (6.7b) 


Ca + 1 


+ 2 + Xq 


77 + 2 
C a + 1 


- 2 + 2 x 3 x|c a + x 2 (xi + x 3 (?7 + 2 )), 


defined on the phase space 

T = {(xi, x 2 , x 3 ) : xi > 0 , x\c a + 2 x 2 + x 3 < 1 } . 


(6.7c) 


( 6 . 8 ) 


The equilibrium points of the system (6.7) are given in table 4. Let us discuss their 
stability. 

1. The equilibrium point P\ is always a saddle. It satisfies K = a? asymptotically, which 
implies ei 1 ~ e~ T = i, e 2 2 ~ e~ T = K Since it <C a as r — > 00, it follows that N <C r~ l . 

2. Although the equilibrium point P 2 can be an attractor for 77 > l,c a < —77 — 3 or 
77 > 1, —77 — 3 < c a < —1, since it can never belong to the phase space (denoted ^ T 2 
in table), we do not discuss it further. 


3. The equilibrium point P 3 is always a saddle. 

4. The equilibrium point P 4 is a source for 77 > 2 , 2^ < c a < ^(77 — 2 ) (77 + 6 ). Otherwise 
it is a saddle. 
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Label 

Xl 

X2 

2?3 

Existence 

Pi 

0 

0 

0 

always 

P 2 

— 2 — r] 

0 

1 

^ T 

P 3 

0 

_ 2 

V 

4 (c a +l) 

T \ 2 

Ca < Ai 

Pa 

0 

V+2 

(c a +l)((?/-2)(?7+6)-16c a ) 


— Aca+IJ—2 

(-4c a +?7-2) 2 

^ a, lA2 ^ 4 





> 2, ^ < C a < A 2 

P 5 

?y(ca+l)( 4 c a + 3 ) 

1 

(c a + l)( 4 c a + 3 ) 

?7> 1,-1 <C a < -| 

(2c a +l) 2 

—2c a —1 

(2c a +l) 2 

Pe,7 

0 

1 

l±-\/l+ c a 

0 

always 


Table 4. Equilibrium points of the system (6.7). We use the notation Ai = - (rf + 4?y — 4) and 
A 2 = ^ (? 7 2 + 477 - 12 ) . 


Label 

Eigenvalues 

Pi 

P 2 

P 3 

Pa 

P 5 

P(i,7 

-1,2,2 

O 2 | Ca+77+3 Ca+^+3 

\Z(Ca“hl) (Ca+77+3) -^(Ca + l) (Ca+77+3) 

0 r/+2+^/ 64c a —7 ?j(?7+4)+36 ri+2—y] 64c a —7ry(r;+4)+36 

277 > 277 

»? 2 -4 o 77 ( 77 + 2 ) 0 77 ( 77 + 2 ) 

2(4ca-r?+2) "b 4 c a - 7 j +2 4c a -?j+2 

V 0 1 0 

2ca+l’ -2co-l 

2((c a = p-y^ c a - l - l)H _ l) ( 77 + 6 ) (^lzh\/ c a +l)+4c a 2^2c a = F\/ca+l)H _ l) 

c <* ’ (Vca+l-l) 2 ’ c « 


Table 5. Eigenvalues of the equilibrium points of the system (6.7). 


5. The equilibrium point P§ is a saddle for rj > 1, —1 < c a < — It is non-hyperbolic for 
c a = — | [but it behaves as a saddle]. 

6. Pq is a source for rj > 1, — 1 < c a < 0, or rj > 1, c a > 0. 

7. Pj is a sink for 1 < 7 / < 2, -P(j] — 2)(rj + 6) < c a < 0. It is a source for 1 < 77 < 2, c a > 0 
or r) > 2, c a > ^(77 — 2) (77 + 6). It is a saddle otherwise. 


6.2 Static vacuum aether with a scalar field with harmonic potential. 

Let us investigate a static vacuum aether with a scalar field with harmonic potential V (</>) = 

2 / 2 

(also see [15]). Introducing the new dimensionless variables 


u K \/2m 

X2 1 0 •> ^5 • Xq 

a a z a 


e i(^ > ) 

a 


x 7 


mcf) 

V2 a 


(6.9) 
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we obtain the dynamical system 


dX 2 = 2x\c a - — — ^y 1 + x 2 (x 4 + 4 x 2 + x\ - 2) , 


= 2X4 2X2 c a — 


dr 
dx 4 
dr 

'y— = X 5 ( 2 x 2 C a - 

dx 6 
dr 
dx 7 


Ca + 


r 2 

2 x 7 i i o , 2 


C a + 1 


+ X 4 + 3X2 4“ ®g — 1 I j 


C a + 1 


+ X4 + 3x 2 + Xq 


= X 6 


C a + 1 


+ x 4 + 4 x 2 + Xg — 2 + 2x2X 6 c a — X5X7, 


dr 


= x 7 


C a + 1 


+ X4 + 3 x 2 4 " Xg I + 2x2X7 c a — 


X 5 X 6 


subject to the constraint 


1 = c a x 2 4- 2x 2 4- x 4 + -Xg - X 7 . 


( 6 . 10 a) 

( 6 . 10 b) 

( 6 . 10 c) 

(e.iod) 

( 6 . 10 e) 

( 6 . 11 ) 


The constraint (6.11) is preserved by the dynamical system (6.10). Solving the constraint 
( 6 . 11 ) for X 4 and substituting back into the system ( 6 . 10 ) we obtain the reduced system: 


dx 2 

dr 


= X 2 C a - — -+ ^X 2 (4x 2 + 2 X 7 + Xg - 2) , 


c a 4- 


dx 5 _ 1 . 

dr ~ 2 Xb ( Ca 


c a 4-1 


4- x 2 ) + 2 x 2 4~ Xg 4-21, 


dx g 
dr 
dx 7 
dr 


X6Ca ( c X \ 1 + X2 ) + \ ( 4x2 + x l ~ 2 ) - 2x 5 x 7 ) 


= x 7 c a 


X 7 


c a 4-1 


+ x 2 — 


X 5 Xg 


Xa 


+ X7 X 2 + - 7 T + 1 


defined in the phase space 


T = \ (x 2 , x 5 , Xg, x 7 ) : c a x 2 + 2 x 2 + x\ 4 - -x@ - X 7 < 1 \ . 


( 6 . 12 a) 

( 6 . 12 b) 

( 6 . 12 c) 

( 6 . 12 d) 

(6.13) 


The equilibrium points of the system (6.12) are described in tables 6 and 7. Let us 
discuss their stability. 

1. Q\ is always a saddle. 

2. The line of equilibrium points Q 2,3 is normally hyperbolic and is stable when X 2 > 2. 

3. The equilibrium points are non-hyperbolic. They have a 3D stable manifold and 
a ID center manifold for — 1 < c a < — ^ and a ID stable manifold and a 3D center 
manifold for c a = — |. 

4. Qaj are non-hyperbolic. They have a 3D stable manifold and a ID center manifold 
for — ^ < c a < —7j. They have a 3D unstable manifold and a ID center manifold for 
c a < — 1 or c a > — ^ [the non zero eigenvalues are always of the same sign]. 
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Label 

X2 

x 5 

x 6 

x 7 

Existence 

Qi 

0 

0 

0 

0 

always 

Q2,3 

x 2 

0 

1 — 2 x 2 ~ CaXlf 

0 

1 — 2X2 ~ °a x 2 2 > 0 

Q 4,5 

Q&j 

1 

0 

0 

0 

, \/2\J (c a +l)(4c a +3) 

1 \/ (-Ca-l)(4c a +3) 

(c a + 1) (4:C a + 3) < 0 

(c a + 1) (4:C a + 3) > 0 

2c a +l 

1 

± 2c a + l 

0 

2Ca + l 

± 2c a +l 

Q 8,9 

1 

0 

0 

0 

always 

l±-\/ C a + 1 


Table 6. Equilibrium points of the system ( 6 . 12 ). x\ is a parameter and hence the curves 
represent lines of equilibrium points {x\ = 0,^2 = 2 are special points on these curves). 


Label 

Eigenvalues 

Qi 

- 1 ,- 1 , 1,1 

Q2,3 

2 -x£, 2 -x$, 2 (l -s 2),0 

Q 4,5 

4c a +3 4c a _ l _ 3 o q 

2ca+l> 2co+l> Z ’ U 

Q 6,7 

4(c a +l) 4c a +3 4c a +3 n 

2c a +l ’ 2c a +l > 2c a + l ’ 

Q 8,9 

(2c a = FV c a - l - l)+l (2c a = F\/ c a + l)H _ l 2^C a =F\/Ca + l ) + l) p 

Ca > Ca ’ Ca ’ U 


Table 7. Eigenvalues of the equilibrium points of the system ( 6 . 12 ) given in the previous table. 

5. Q 8 is non-hyperbolic. It has a 3D unstable manifold and a ID center manifold for 
— 1 < c a < 0 or c a > 0. Otherwise, its center manifold has dimension greater than 1 . 

6 . Qg is non-hyperbolic. It has a 3D stable manifold and a ID center manifold for —j< 
c a < 0. It has a 3D unstable manifold and a ID center manifold for c a > 0. Finally, 
Q 9 has a 2D unstable manifold, a ID center manifold and a ID stable manifold for 
-1 < c a < -§. 

Although static models are of particular physical importance, in this paper we have 
primarily focused on the mathematical properties of the solution space. It can be observed 
that the phase spaces ( 6 . 8 ) and (6.13) are in general non-compact; thus a more detailed 
analysis requires the introduction of compact variables. In addition, since the equilibrium 
points in (6.13) are generically non-hyperbolic, the use of the center manifold theorem is 
required, which is beyond the linear analysis provided here. A more detailed stability analysis 
for the equilibrium points of both the dynamical systems (6.7) and (6.12), and the study of 
the tilted aether static model, is left for the companion paper [46]. 
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7 Discussion 


In this paper we have studied spherically symmetric Einstein-aether models with tilting 
perfect fluid matter, which are also solutions of the IR limit of Horava gravity [ 8 ]. We used 
the 1+3 frame formalism [19-21] to write down the evolution equations for non-comoving 
perfect fluid spherically symmetric models and showed they form a well-posed system of 
first order PDEs in two variables. We adopted the so-called comoving aether gauge (which 
implies a preferred foliation, the only remaining freedom is the coordinate time and space 
reparameterization freedom). We also introduced (/?-) normalized variables. The formalism 
is particularly well-suited for numerical and qualitative analysis. In particular, we considered 
the special subset U = v = 0 (where we also assumed cq = 0 and c a / 0) and derived the 
final reduced phase space equations in normalized variables. 

The formalism adopted here is appropriate for the study of the qualitative properties 
of astrophysical and cosmological models with values for the non-GR parameters which are 
consistent with current constraints. In particular, motivated by current cosmological obser¬ 
vations, we have studied inhomogeneous cosmologies in Einstein-aether theories of gravity. 

We first considered dust models. We investigated a special dust model with U = 0 and 
v = 0 in normalized variables (assuming c a / 0 ) and derived a reduced (closed) evolution 
system. The FLRW models in this special dust model correspond to an equilibrium point. In 
these models we are particularly interested physically in their late time evolution. Therefore, 
we paid particular attention to the sinks for different values of the parameter c a (which were 
summarized earlier). In all cases O, —>• 0 to the future. For all solutions with small c a < 3/8, 
T> —» 1 (Q —>• 1) and the shear goes to zero at late times. Consequently, the models close to 
GR isotropize to the future. 

We briefly reviewed the FLRW models in which the source must be of the form of a 
comoving perfect fluid (or vacuum) and the aether must be comoving. We then considered the 
spatially homogeneous Kantowski-Sachs models [36] using appropriate normalized variables 
(which are bounded; note that we did not use the /3 —normalization here), and obtained the 
general evolution equations. We then considered a special case with 3cq = C 1 + 3 C 2 +C 3 = 0 and 
analysed the qualitative behaviour. In this case a full global dynamical analysis is possible, 
and we determined both the early and late time behaviour of the models and their physical 
properties. 

In the case c a < \ (i.e., c > 0), when 7 < |, there is the unique shear-free, zero 
curvature (FLRW) inflationary future attractor (P2), and for | < c a < \ (i.e., 0 < c < ^) 
and 0 < 7 < 2 all of the sources and sinks (respectively, P4 & Ps and P 5 Sz P?) have maximal 
shearing and all except one sink (T 7 ) have zero curvature. For c a < | (i.e., c > |), the 
points TV & P% do not exist, and the sources and sinks with maximal shearing are P 4 and 
T 5 , respectively. 

Finally, we considered static models for a mixture of a (necessarily non-tilted with v = 0) 
perfect fluid with a barotropic equations of state and a scalar field (with a self-interaction 
potential V that depends only on the scalar field). In particular, we studied the special cases 
of a tilted perfect fluid and no scalar field (previous work had assumed a comoving aether and 
a comoving fluid) and a stationary vacuum with a scalar field (with a harmonic potential). 
The equilibrium points in the resulting dynamical systems in these two cases were determined 
and their stability was investigated. Although models that are not evolving with time are 
of physical importance, and physical information can be obtained from their qualitative 
analysis, we have primarily focussed on the mathematical properties of the solution space in 
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this paper. The physical interpretation of this analysis will be comprehensively discussed in 

[ 46 ]- 

We also examined the conditions for the existence of McVittie-like solutions in the 
context of Einstein-aether theory. We found that they only exist for the choice of parameters 
c a = 0,7 = 0, and for an aligned aether (v = 0). Since 7 = 0, the matter fluid corresponds 
to a cosmological constant (and 6 is always a constant). Irrespective of the sign of the initial 
expansion, the physical variables tend to zero as t —>• +oo. 

In future work we shall investigate the general Kantowski-Sachs models and the static 
models more comprehensively. In particular, it would be of interest to determine the structure 
of stationary rotating solutions; rapidly rotating black holes, unlike the non-rotating ones, 
might turn out to be very different from the Kerr metrics of GR. 

We note that the tilt is defined relative to matter; one important question is to in¬ 
vestigate whether this tilt decays to the future in general. We shall also study spherically 
symmetric, self-similar spacetimes which also admit, in addition to the three Killing vectors, 
a homothetic vector [ 21 ]. 
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A Models and the parameters c* 

We can study different models with different dimensionless parameters Cj. From earlier: 

eg = C2 + (ci + c 3 )/ 3 , c a = ci + C3, Cuj = ci — C3, c a = C4 - ci. 

Since the spherically symmetric models are hypersurface orthogonal the aether field has 
vanishing twist and the field equations are therefore independent of the twist parameter c u 
[ 8 ] (this is equivalent to being able to set C 4 = 0 without loss of generality [ 1 ]). 

A second condition on the Cj can effectively be specified by a renormalization the New¬ 
tonian gravitational constant G. From [1] we have that G/v = G (l — ^(ci + C4)) . So long 

as (ci + C4) < 2, so that the gravitational constant is positive, we can effectively renormalize 
and specify (ci + C4). If not, and we reduce the theory to a one parameter model, the theory 
might be pure GR in disguise [in GR c. L = 0]. 

The remaining two non-trivial constant parameters in the model must satisfy additional 
constraints (it will be useful here to define c 2 = 1 — 2 c a < 1 ): 

Observations: The models (i.e., the values of the Cj) must be consistent with all observa¬ 
tions. In general, if the magnitudes of all of the Cj are (non-zero and) small (e.g., less than 
ICR 2 ), then the models will be physical [1, 24]. 

Self-consistency: There are also a number of self-consistency requirements [1, 24]: 


0 < ci + c 3 < 1 , 


(A.l) 


0 < ci 


< (g1 + c 3) 

3 “ 3[1 - (ci + c 3 )] ’ 


(A.2) 


which can be written in terms of eg , c CT , c a . Note that this imples that 0 < ci, C3 < ci < 1— C3. 


A.l Case A: 

All of the Cj are small (and not all zero). We set C 4 = 0. We renormalize and chose ci so 
that ci + 3 C 2 + C 3 = 0 (i.e., eg = 0). The self-consistency relations (A.l, A. 2 ) then imply that 
0 < C 3 < ci < 2 c 3 . We thus have a two parameter model with small c a >0, c a < 0. Using 
ci = dc 3 (1 < d < 2), we have that c a = |(1 — c 2 ) and c a = — (i+^ Co-. Note that it is not 
possible for c a = —c a in this case. 

Summary case A: c a = ^(1 — c 2 ) > 0, c a = — c c < 0, eg = 0 

A.2 Case B: 


We set 


c 4 = -- 


C3 

ci 


c 2 = -- 


( 2 ci 2 + cic 3 - c 3 2 ) 
3ci 


(A.3) 

(A.4) 


(before the field redefinition of C4), so that the parameterized post-Newtonian (PPN) param¬ 
eters ai = 0:2 = 0 (and hence all solar system tests are trivially satihed [1]). A two parameter 
family of models (ci 7 ^ 0 , 03 ) satisfying the consistency conditions (A.l, A.2) results. Here 

(Summary): 
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and the c* satisfy 


^>‘ 2 + cs2) <0. 

(A.5) 

Cl 

0 < Cfj = Cl + c 3 < 1, 

(A.6) 

(Cl 2 - C 3 2 ) ^ 

ce = „ <0, 

3ci 

(A.7) 

— 9 Cq 2 ) = (c a + 3 Cq + 2 Co-) . 

(A.8) 


A.2.1 Case B(ii): 

If we also renormalise the Newtonian gravitational potential by setting ci + C 4 = 0, then 
from (A.3) we obtain ci 2 = C 3 2 , and hence ci = C 3 (since ci = —C 3 cannot satisfy (A.l, 
A. 2)). In this case cq = 0. Hence, c a = — c a , and we have that c a = ^(1 — c 2 ) > 0 and 
c a = — 4 ( 1 “ c 2 ) < 0 , where c 2 < 1 but need not be small (and the self-consistency relations 
(A.l, A.2) are satisfied). 

Summary case B(ii): c a = |(1 — c 2 ) > 0, c a = — ^(1 — c 2 ),cq = 0 

A.3 Case C: 

In principle we can study the physics of the models for different parameter ranges of the C{. 
If we study the models in the early universe (where the constants c* can be replaced with 
evolving parameters [5]), then the observational constraints above need not apply. 

In one particularly interesting theoretical case (see section 4), we could consider 3cg + 
2 c a = 0 [we could also use the renormalization of the Newtonian gravitational constant and 
consider the case c a = 0]. Note that 3c# + 2c CT = ci + C 2 + C 3 = 0. This implies that 
the PPN parameter a\ diverges [1] and conditions (A.l, A. 2) can only be satisfied when 
ci = C 2 = C 3 = 0, the GR case. However, for theoretical reasons it may be of interest to 
study this case in early universe cosmological models. 

Summary case C: c a = ^(1 — c 2 ) > 0, eg = —1(1 — c 2 ) < 0, [ c a = 0] 
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B Further development of the governing equations 


Let us further develop the governing equations presented in Section 3. Combining equations 
(3.13d) and (3.13e) and using the identity £ = 1 — Q, we obtain 


3/C 1 ATT 

~ 2(2 c„ - 1) + 3 A 


+ ^r(c a + 1) 


1 


+ 


2 (c a + 1) 1 — 2c a 

3cg + 1 + 

1 


1 — 2c„ 


+ 


A 2 


2 ( 2 c a - 1) 


3 cq + 1 2c a — 1 

4 (3 c e + 2c g ) Q 

3cg + 1 


1 - 2c, 


(c a + 1 ) ( 3 Cq + 2c a )d\U 
3 (3 c e + 1) (2c a - 1) 

+ 


+ 


JJ }_jj 2 ( 5c a + 2 2 (c a + 1) 


3 Cq + 1 


6 — 3c0 — 14 c ct + 2 - 4 c a + ^ 


4c a -2 


3 Cq + 1 


Cl (u 2 (3(7 + l)cfl + 2(7 - 2)c a + 3 ) — 67 c a - 3c<9 + 4 cq + 37-3) 


2 (3c0 + 1) (2co- — 1) (1 — v 2 ) 


+ 


+ 


1 — 9 Cq — 8 c a 
6 ce + 2 


(B.l) 


Combining equations (3.13c), (3.13f) and (3.14d) and (3.14e) and using the identity £ = 1 —Q, 
we obtain 


r = (q- l)c a U - c a d 0 U - 3A(1 - Q) + 


3y vQ 
2(1 - v 2 )' 


d 0 A =(q + 3(1 - Q))A -U + r, 


diQ = - 


3 (2 Cg - — 1) > 4(1 — Q) 

3cg + 2Ccr 


+ r[Q + 


1-2 c„ 


3^yvQ 


3c e + 2c a J 2 (3ce + 2c a ) (1 - v 2 )' 


(B.2a) 

(B.2b) 

(B.2c) 


Finally, solving the equations (B.l) and (B.2a) for q and r we obtain: 


Q = 


1 



2c a (cg+l)(3c0+2c CT )U 2 \ 
3(3cg + l)(2c CT —1) ) 


3^4 2 

1 — 2c ff 


4 (Cg + l) LC4 
3c e + l 


2(1-20 UA _ Rn2 2 (cq + 1) U 2 (5cq + 2) U 2 , 

2c a - 1 2 3 c e + 1 2 c ct - 1 

12(Q — l) 2 (2cq — 1) 3Q(-16cq + Q(3c 0 + 14cq-6) + 8) 

3c e + 1 2 C q - 1 

3 ((7 — 2)u 2 — 3y + 2) fl 3 ((7 + l)u 2 — l) fl 
(3cg + 1) (u 2 - 1) (2c ff - 1) (u 2 - 1) + 

| (Cq + 1) (3c e + 2 Cq) U (~6A(Q - 1) + 2c a (d 0 U + U) + fg) j 
+ (3 c e + 1) (2 C q - 1) + 

2d\U ( C q + 1) 3 9/C 20!LT (c a + 1) 3 

3Cfl + 1 + 1 - 2 C q + 2 C q - 1 + 2 C q - 1 ^ 2 C q - 1 + 


(B.3a) 
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r = 3QA — 3 A — c a doU — c a U + 
c a U 


{■ 


3 r yv£l 

2(1 - v 2 ) + 

3 A 2 4 (ca + 1) UA 


(n _ 2c a (c a +l)(3cg+2c a -)[/ 2 \ l 1 — 2c, 
\ 3(3c 9 +1)(2c ct -1) ) 


3 cq + 1 


+ 


+ 2 ( 1-2 c a )UA _ 6q2 + 2 (c a + 1) U 2 + (5c a + 2 ) U 2 | 


2c„ - 1 


3cg + 1 


2c, - 1 


+ 


12(Q — l ) 2 (2c, — 1) , 3S(-16c <T + Q(3c e + 14c a -6) + 8 ) 


+ 


3 Cq + 1 2 C a — 1 

3 ((7 - 2)v 2 - 37 + 2) n 3 ((7 + l)v 2 - 1 ) n 


+ 


( 3 c e + 1 ) (v 2 - 1 ) 


( 2 c, - 1 ) ( v 2 - 1 ) 


+ 


+ 


S'yvQ \ 


(ca + 1) (3 c e + 2O U (~6A(Q - 1) + 2 c a (d 0 U + U) + 


2d\U {c a + 1) 
3c e + l 


(3c e + 1) ( 2 c ct - 1) 
3 9/C 


+ 


1 - 2c, 


+ 


2CW - 1 


2d\U (c a + 1) 
2 c,- 1 


2c,t - 


i+4 


(B.3b) 


These equations give the expressions for q and r in terms of the normalized variables and 
the derivatives d()U,diU. The expression (B.2c), is used to eliminate the spatial derivative 
d\ Q from the equations. 

The final equations for the reduced phase space ^ 'N,U,E\,KL , Q, are equations 

(3.13a,3.13b,3.13d,3.13g,3.13h) and (3.14e) for do A, subject to the constraints (3.14a-3.14d) 
and a constraint for d\ Q (rather than 5iS), with q and r defined as in (B.3). From the 
equations we have either the evolution equation for U given by (3.13c) or the definition of r 
given by (B.2a) (or by (B.3)(b), after ^-elimination in (B.2a)). 

The commutator equation (2.30) can be expressed in terms of the normalized variables 
by 

[0o,3i] = ^2 [ e o,ei] + (1 + q)d 1 - rd 0 

= (U-r)do + (q + 3-3Q)d 1 . (B.4) 


On the other hand 


d 0 r = d 0 di(-ln/3) = -[3 o ,0i]ln^ - 0i0 o hi/3 
= U(1 + q) + r(2 - 3Q) +d\q, (B.5) 

where we have used 9oln/3 = —(1 + q) and q is given by (B.l). 10 

10 Note that the evolution equation for r will contain the second order derivative term didiU (via the term 
diq), as occurs in the GR setting [ 42 ]. 
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The final equations for the reduced phase space ( AT, U, E\,JC, Q, Q,v,r) are then 


doE\ = (q + 3E)£j ’, 
d 0 lC = 2 qJC, 

c a d 0 U = (q - l)c a U - r - 3.4(1 - Q) + 


2(1 — v 2 )' 


a„Q - ( 7l 1)a ;( > = e( i +„ - a) - r -M + ^ - 2 l c ° +1M , i> + 


3(3« + l) 


3(3c 0 + l) 3(3c e + l) 


(c a + 1) U 2 2 (2cp- — 1) (1 — Q) 2 12(( 7 -2)u 2 -3 7 + 2) 

+ 3 (3c 0 + 1) + 3c 0 + 1 2 (3cg + 1) (1 - v 2 ) ’ 

doA = (q + 3£).4 — U + r, 

Q 0 n + (7 ~ 2 )cffil 2 _ 7^1 v _ 27 AvQ + 2 qQ (1 - (7 - l)u 2 ) | 




G_ 


G_ 


G _ 


3 7 Q (l — u 2 ) 12 ^ 2(7 — 2 )rul 2 212 ((1 — 2 7 )u 2 + l) 


G_ 






(7 — 1 ) (l — u 2 ) 2 di!2 (7 — 2 )vd\v _ 2(7 - 1).4 (l - v 2 ) 


d 0 v - 


7 12G_ 


+ 


G_ 


G_ 


-+ 


3(7 — 2 )Q (l — u 2 ) u 2(7 — l)r (l — u 2 ) 2 2 (l— u 2 ) v . 2 

n 7 G_ + G_ + 1 w ' ’ 


G- 

dor = U(l + q) + r(2- 3 Q) + dig, 


(B. 6 a) 

(B. 6 b) 

(B. 6 c) 


(B. 6 d) 

(B. 6 e) 

(B. 6 f) 


(B. 6 g) 

(B. 6 h) 


where 


3/C 1 ATT 

9 ~ 2(2Co- - 1) + 3 AU 


+ ^r(c a + l) 


1 


+ 


2 (c a + 1) 1 — 2c q 

3c0 + 1 + " 

1 


1 — 2 c„ 


+ 


4l 2 


2 ( 2 c ff - 1 ) 


3 eg + 1 2c a — 1 

4 (3cg + 2c ct ) Q 

3cg + 1 


1 - 2 c, 


(c a + 1) (3 eg + 2c a )d\U 
3 (3c e + 1) (2c ct - 1) 

+ 


+ 


JJ 1 fj 2 5c a + 2 2 (c a + 1 ) 


3cg + 1 


6 — 3cg — 14c ct 2 — 4c ct \ _ 9 
+ -- 7~T + 1 Q + 


4c ct — 2 


3cg + 1 


12 (u 2 ( 3(7 + l)c e + 2(7 - 2)c CT + 3) - 67 C 0 - 3 Cg + 4co + 3 7 - 3) 


2 (3cg + 1 ) ( 2 c ct — 1 ) (1 — v 2 ) 


+ 


+ 


1 — 9 eg — 8 c 0 


6 cg + 2 
subject to the constraints 

d\M~ l = {r-U)N~ l , 
diK. = 2 (r + A)1C, 

d\A + c a d x tj = ——AC + r (^c a U + .4^ + 2c a AU — —c a U 2 + —.4^+ 

- | (3c0 + 2co) Q 2 + 3 (2co - 1) Q + ^ ^ ~ l)V l + ^ - 3co + 

2 2 (1 — ir) 2 


aic = — 3 ( 2 %- 1 |f- Q) +.-|g + 

3c^ + 2 c a 


2(1 — u 2 ) 

1 — 2c n 


37^12 


3 c 0 + 2 co/ 2 ( 3 c 0 + 2 co) (1 - v 2 )' 


(B.7) 


(B.8a) 

(B.8b) 


(B.8c) 

(B.8d) 
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We can choose to study the system (3.13), (3.14), (B.3) or the system (B.6), (B.7), 
(B.8), depending on the particular application. 

In the comoving aether temporal gauge, which implies a preferred foliation, the only 
remaining freedom is the coordinate rescalings t —>■ f(t ) and x —>• g(x) (time and space 
reparameterization freedom), consistent with 

d 0 El = (q + 3X)El, 
diAT 1 = (r-U) AT 1 , 

where we recall that do := AT~ 1 dt, d\ := E\d x . 

How do we best treat in the evolutions eqns? There is no evolution eqn for Af~ l 
but, in principle, we can integrate the spatial constraint and use the time reparameterization 
to determine A A -1 . [This is what happens in GR in the separable gauge, where neither an 
algebraic equation or an evolution equation for Af" 1 is available; the conditions are sufficient 
to integrate the spatial derivative constraint equation, whence a time redefinition can be 
employed to set AT~ 1 = 1.] 

Depending on the application, we could do one of the following: (i) Since A is positive- 
definite, we can determine the qualitative behaviour of the system by simply studying the 
right-hand-sides of the evolution equations, (ii) In many special cases of interest we can 
integrate for Af^ 1 and replace the left-hand-sides by partial time derivatives (as in GR in 
the separable gauge), (iii) Numerically we could solve for AT -1 in the integration (although 
this looks messy), (iv) Analytically, we could, for example, use the commutators to obtain 
evolution and constraint equations for r, and then Af^ 1 = exp (f (U — r)dx) and define new 
variables (e.g., use U — r as a variable), (v) While the comoving aether gauge choice is 
motivated physically, it is not ideal for doing analysis and numerics; we could change to a 
separable gauge. 

In practice, it is often useful to choose a gauge in order to compare with the FLRW 
model as easily as possible (e.g., so we can choose integration functions for A f,E\ in the 
above to be trivial). 

As an example, let us consider McVittie-like models [43]. The line element is given by 
) dt 2 + A 2 ( 1 + 1 [dx 2 + x 2 (dd 2 + sin 2 0d<p 2 )\ , (B.ll) 

1 + T47/ V 2Ax J 

where A = A(t ) and M is a constant. We recall that as x —» +oo the metric approaches the 
flat FLRW solution and for constant A we obtain the Schwarzschild solution. Several aspects 
of the McVittie solution, of geometrical and physical relevance have been investigated by 
many authors [44], In GR the McVittie solution is unique under the following assumptions: 
(i) The metric is spherically symmetric with a singularity at the centre, (ii) The matter 
distribution is a perfect fluid, (iii) The metric must asymptotically tend to an isotropic 
cosmological form, (iv) The fluid flow is shear-free. McVittie-like models were investigated, 
e.g., in the references [45]. 



(B.9) 

(B.10) 
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For this metric we obtain in our scenario: 


which implies that 


16Mx 2 A 2 

7/ = - 

(M — 2 xA)(M + 2 xA) 3 ’ 
AxA(M — 2xA) 
a ~ (M + 2xA ) 3 ’ 

_ 16 x 2 A 2 

X ~ ( 2 xA + M) 4 ’ 

3 i? = 0, 

3 64 Mx 3 A 3 

S+ = ( 2 xA + M) 6 ’ 


£j+ = 0, 


0 = 


3A 

T’ 


(B.12a) 

(B.12b) 

(B.12c) 

(B.12d) 

(B.12e) 

(B.12f) 

(B.12g) 


Q = l,E = 0,r = 0, 

_ AA( 2 xA + M) 2M 
Q ~ A 2 {M — 2xA) M-2xA 


For this metric, the equations (B.6c) and (B.8d) lead to a contradiction, unless c a = 0 and 
either v = 0 or 12 = 0, or both. Assuming c a = v = = 0 and substituting into equation 

(B.6d) we obtain 

AA(2xA + M) — 2MA 2 = Q 

A(M- 2 xA) 1 

which, in general, is not satisfied as we vary t and x. Thus, assuming c a = 0 and v = 0, ^ 0, 
we obtain 

2 (3 c e + 1) (AA( 2 xA + M) - 2 MA 2 ] 

n = -4- ; - L, (B.14) 

(3y — 2)A 2 {M — 2xA) 

The equations (B.6) and the constraints (B.8), with the exception of (B.6f), (B.6g) and 
(B.8c), are identically satisfied. The equations (B.6f) and (B.6g) reduce to 


d 0 Q = ( 2 q + 2 - 37 ) 11 , 
( 1 - 7 ) 01 ^ = 7 ^ 17 . 

which, after the substitution of (B.14), lead to 
A 2 (( 7 + 1 )M — 2(7 — 1 )xA) 


A = 


A = 


A( 2 xA + M) 

A 3 (4 xA ((7 - 1)(3 7 - 1 )xA - 3 7 2 Af + M) + ( 7(37 + 2 ) + 1 )M 2 ) 


A 2 (2xA + M ) 2 

and the restriction (B. 8 c) becomes 

3 (3 c g + 1) (A 2 ((4 - 67 )xA + (37 + 2 )M) - 2 AA( 2 xA + Af)) 
2(37 — 2)A 2 {M — 2xA) 


(B.15a) 

(B.15b) 


(B.16a) 

(B.16b) 


= 0 . 


(B.17) 
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The last three equations, in general, are not satisfied simultaneously for all t and x unless we 
set 7 = 0 , which implies 

A(t) = C2e cit , 0 = 1 + 3c0. (B.18) 

Summarizing: the McVittie-like models only exist for the choice of parameters c a = 
0,7 = 0, and for aligned aether (v = 0). Since 7 = 0, the matter fluid corresponds to a 
cosmological constant. The solution is characterized by 

16cnMx 2 e 2cit 

^ —-—- 

(M — 2c2xe cit ) ( 2c2xe cit + M) 3 ’ 

4c2xe Clt (M — 2c2xe Clt ) 

( 2c2xe Clt + M ) 3 

16 c\x 2 e 2cit 

(! 2c2xe Clt + M) 4 ’ 

3 R = 0, 

3 64:C 3 Mx 3 e 3cit 

+ ( 2c2xe cit + M) 6 ’ 

o+ = 0, 

9 = 3ci. 

and 

Q=l,S = 0,r = 0,u = 0,O = l + 3 c e , q = -1. (B.20) 

Irrespective of the sign of the expansion (9 > 0 or 9 < 0), the quantities u,a, K, 3 S + 
tend to zero as t —>• +oo. 


(B.19a) 

(B.19b) 

(B.19c) 

(B.19d) 

(B.19e) 

(B.19f) 

(B.19g) 
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B.l Special case: u = 0 

The evolution equations above were derived under the assumption that a / 0 
consider the special case of u = 0 and display the appropriate eqns: 


eoCei 1 ) = -\{0 - 6cr + )ei 1 , 
e 0 (K) = ~l(9 + 3a + )K, 

_ (p , 1 P 2 . 6(2^-1)4 (2-3 7 +( 7 -2 W) 9 

eo(9) - “ + + 3c+ 1 + 2(3c + 1)(1 * « 2 ) ’ 

, , 1 2 . {3c e + l)9 2 a 2 K 

+J 2 + + 18(200- — 1) 2 ( 2 c ct - 1) 2 ( 2 c ct - 1) 


+ 


((1 + 7)u 2 - l) jx 

6 ( 2 c<j — 1)(1 — v 2 ) ’ 


e o(a) = ~^ a (0 + 3 cr+) + 

_ e i (A) v(2-'y) _ 7 /rei (v) = 2^{xv 2 a + + 7 {v 2 ~ 3) A# _ 


G _ 


G_ 


G_ 


3G- 


G _ ’ 


ei (/i) (l - u 2 ) ( 7 -I) u (2 - 7 ) ei (u) 2 u (l - v 2 ) a + 
eo[v) -—-^-=-—- 1 - 


7 9 G- 


G _ 


G_ 


v (l — v 2 ) (37 — 4) 9 + 2v 2 (l — v 2 ) (7 — 1) a 


3G- 

Constraints (3 eg + 2c a / 0): 




ej (In K) = 2a, 

/ \ (j_j -LL 1 , . 9 3 . .9 lii 3a 

6l(a) = 2(1 — v 2 ) ~ 6 ( 6 + )0 ~2 ( " 1} a + ~ ~2 + ~ 2 ~ 


eiW = - 


37 c a vfi 


ei(cr + ) = 3 acr + — 


(3c 0 + 2c< r )(l - v 2 ) ’ 
3^cevfi 


2(3cq + 2c ff )(l — v 2 )' 


We now 

(B.21a) 

(B.21b) 

(B.21c) 

(B.21d) 

(B.21e) 

(B.21f) 

(B. 21 g) 

(B.22a) 

(B.22b) 

(B.22c) 

(B.22d) 


- 45 - 






















B.1.1 Normalized variables 

The /3-normalized equations are: 


d 0 E\ = {q + 3Z)E\, 
d 0 IC = 2 qJC, 


(B.23a) 

(B.23b) 

(B.23c) 


1 „ VC A 2 

= S(1 + q — 3 Q) + -£ + —— — — — — -r —jz tt 

2 2(2 c CT -l) 2(2ccr — 1) 2(2co- — 1) 


, (3eg + 1) Q 2 | 


+ 


fl ((7 + l)u 2 - l) 

2 ( 2 c ct — 1 ) (1 — v 2 ) 


doA — qA T 


3q vCl 


doCl + 


2(1 — v 2 ) ’ 

(7 — 2 )vd\Cl ^Cldiv _ 2 ^AvCl 2 qCl (l - (7 - l)u 2 ) 


(B.23d) 

(B.23e) 


G- 




G_ 


G_ 


+ 


3 7 g(l-u 2 )ft 2( 7 -2)n;fi 211 ((1 - 2 7 )u 2 + l) 

+ G~ + G~ ’ 

(7 - 1) (l - v 2 ) 2 did ('y — 2 )vdiv 2(7 — l)A (l - v 2 ) 

doV - ^gZ - + - GZ - = - GZ - 

3(7 — 2)Q (l — v 2 ) v 2(7 — l)r (l — v 2 y | 2 (l — v 2 ) v 

gZ ^gZ + gZ ’ 

subject to the restrictions: 

d l A/" 1 = rU~\ 
diK = 2(r + A) 1C, 

d\A = — —1C + rA + — A 2 — — (3 eg + 2co-) Q 2 + 3 (2ccr — 1) Q 

, 3 H ((7- l)u 2 + 1) o 3 
+ 2(1— v 2 ) ff+ 2’ 


(B.23f) 


-+ 


diQ = rQ - 


37 c a vCl 


(3c e + 2c a ) (1 - u 2 )’ 


d\T> = rS + 3.4.S — 


97 cgvCl 


2 (3 c e + 2 ca) (1 - v 2 )' 


(B.23g) 


(B.24a) 

(B.24b) 


(B.24c) 

(B.24d) 

(B.24e) 


Using the identity S = 1 — Q, and combining equations (B.23c) and (B.23d) and equa¬ 
tions (B.24d) and (B.24e), respectively, we obtain: 


Q = ~ 


3/C 


■ + 


-4 2 


2 ( 2 c ct - 1 ) 2 ( 2 c a - 1 ) 


(3cg + 2c g ) (3 c g + 8 c a - 3) Q 2 4 (3 c e + 2c CT ) Q 


■+ 


2 (3c0 + 1) (2C(j — 1) ' 3c e + l 

Cl (y 2 ( 3(7 + l)c 0 + 2(7 — 2) Co- + 3) — 67C0 - 3 c 0 + 4co + 37 - 3 ) 1-9 eg — 8 c, 

2 (3cg + 1) (2co - 1) (1 - v 2 ) + 


2(3c 0 + l) ’ 

(B.25) 


11 Strictly speaking, we assume v ^ 0 here, since in the dust case below v = 0 leads to a contradiction when 
u = 0 . We shall study an exceptional case later. 
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(B.26) 


The final equations for the reduced phase space (jV, E] .1C. Q,£l,v) T are the evolution 
equations and restrictions displayed above (less the equations (B.23d),(B.24e) for the frame 
derivatives of £), where q and r are defined by (B.25) and (B.26), respectively. (It will be 
useful to define T> = A 2 — 3/C in some computations). Again we note that formally the same 
equations can be obtained from the previous case (ft / 0) by setting u = 0; but here we have 


derived the equations properly. 

B.2 The subset U = v = 0 

Let us consider the special subset U = v = 0. We also assume cq = 0 and c„ ^ 0. We first 
assume that 7 7 ^ 1. The final equations for the reduced phase space are then: 

d 0 El = (q + 3(l-Q))El, (B.27a) 

d 0 IC = 2qJC, (B.27b) 

do Q = Q(l + q-Q) + 2 (2 Co- - 1) (1 - Q ) 2 + hi (2 - 37 ), (B.27c) 

d 0 A = qA, (B.27d) 

d 0 n = (2 q - 37 Q + 2)n, (B.27e) 

subject to the restrictions: 

d\Af _1 = rAf~ 1 , (B.28a) 

Si/C = 2(r + A)K, (B.28b) 

diA = -\K + rA + ^A 2 - 3c a Q 2 + | (2c a - 1) ( 2 Q - 1 ) + (B.28c) 

diQ = rQ, (B.28d) 

diQ = 2 rQ, (B.28e) 

where q and r are defined by: 


q = 1 -3IC + A 2 - 2 c a ( 8 c a - 3) Q 2 + 16c a (2 Cfr - 1 )Q 

2 ( 2 c a - 1 ) I 

- n (— 20,(37 - 2) + 3(7 - 1)) + (1 - 80(2^ - 1)}. (B.29) 

r = -3A(1 - Q). (B.30) 

We recall that doE\ = (g + 3(l — Q))Ej, d\M _1 = rAf~ l , where do ■= M~ l dt , d\ := E\d x . 
[The only remaining freedom is the coordinate rescalings t —>• f(t) and x —> g(x)]. 12 

B.2.1 The special case 7=1 

In the special case of dust, we lose the equation d\fl = 2 rfl, but all of the remaining equations 
are valid with 7 = 1 . 

12 A 2 and K, only appear in the evolution equations via the combination T> = A 2 — 3 /C (but also appear in 
the constraints; e.g., via r ). 
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C Lemaitre-Tolman-Bondi model 


The Lemaitre-Tolman-Bondi (LTB) model in GR [49] is the spherically symmetric dust so¬ 
lution of the Einstein equations which can be regarded as a generalization of the FLRW 
universe. LTB metrics with dust source and a comoving and geodesic 4-velocity constitute a 
well known class of exact solutions of Einstein’s field equations [50] . 

The line element for a spherically symmetric comoving dust is: 

ds 2 = —dt 2 + [-^ dx 2 + R(t , x) 2 [d0 2 + sin 2 ddcj) 2 ] (C.l) 

1 + 2 E 

(where an overdot denotes a f-derivative and a prime denotes an x-derivative). The geometric 
variables, {H, a + , 3 R , 3 5 + }, are defined by 


H + a + 


R 

R 1 


H -2a+ 


3 P 4 (ER)\ 3 , 
R r R? R' ’ 1 


_ J_3 R _ 

12 R? ’ 


(C.2) 


where 

3 R = 4eia — 6a 2 + 2 K, 3 S + = —-ei a + -K, 

O O 

= VTTE _i_ 

R(t, x) ’ R(t, x) 2 


(C.3) 

(C.4) 


The matter variable is defined by p = The 3 free functions R(t,r), E, M are 

subject to the Gauss constraint. We can introduce Hubble-normalized variables [51]: 




n = -?- 
H ’ 3H 2 ’ 




J R_ o 3 s+ 

6 H 2 ’ + 3H 2 ’ 


(C.5) 


where (the deceleration parameter) q = 2T?\_ + and D = — 65+. By introducing a new 

time variable r(t,r) defined by dt/dr = l/H, where H > 0, H decouples from the evolution 
equations. Using these variables, expressed in the coordinates above, the evolution equations 
for GR are (the dynamical system) (see also [52]): 13 


d r Q = (H + 4E 2 - 1)J1, (C.6a) 

<9 r E+ = (in + 2 T? + - 2)E+ - i(l - n - E \-D), (C.6b) 

d r D = (n + 4E^_ - 2E + )D. (C.6c) 

The flat FLRW equilibrium point is given by (f l,E +1 D) = (1,0,0), where (Qfe,5+) = (0,0). 
It is a saddle with eigenvalues (1, — |, 1) [51]. 


13 


We note that the exact GR LTB is not a solution to the Einstein-aether equations (see section 4 ). 
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